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1.  .INTRODUCTION 

The  purpose  of  this  report  is  to  document  the  computer 
programs  that  have  been  written  to  implement  the  prediction 
algorithm  described  in  BBN  Report  3 6‘u3  by  M.  Moll,  R.  M.  Zeskind 
and  W.  L.  Scott  entitled  "An  Algorithm  for  Beam  Noise  Prediction". 
The  authors  assume  that  the  reader  is  already  familiar  with  the 
contents  of  that  report.^' 

1.1  Ambient  Noise  Prediction  Computer  Programs 

The  algorithms  for  calculating  the  probability  density  of 
the  averaged  noise  power  at  the  output  of  a beamformer  in  a 
specified  passband  resulting  from  ship  traffic  in  an  acoustic 
basin  have  been  programed  in  FORTRAN  for  a CDC  660  digital 
computer.  The  algorithms  are  implemented  as  three  separate 
programs  and  their  associated  data  files.  Figure  1 presents 
a block  diagram  of  the  data  flow  of  the  ambient  noise  predic- 
tion computer  implementation . 

The  input  information  is  divided  into  three  types. 
Transmission  loss  data  at  a given  sensor  location  is  read 
into  the  computer  and  stored  on  a file  as  a table  of  trans- 
mission loss  as  a function  of  range  and  bearing.  One  such 
transmission  loss  file  is  created  for  each  frequency  of 
interest.  Sensor,  route  and  ship  traffic  information,  which 
is  independent  of  frequency,  is  read  into  the  computer.  This 
frequency-independent  information  is  used  as  the  input  to  a 
program  which  computes  such  quantities  as  route  segment  lengths, 
earth-centered  angles  between  routes  and  sensor,  and  other  geo- 
metric parameters  needed  in  the  computation  of  the  character- 
istic function  <fy(a)).  The  outputs  of  this  geometry  program 
are  stored  on  a file  for  later  use  as  inputs  to  the 
characteristic  function  program. 

The  characteristic  function  of  beam  noise  power  is  com- 
puted via  the  algorithm  derived  in  Section  2.0  of  BBN  Report 
3653.  This  program  reads  in  radiated  noise  data  for  each 
type  of  ship  and  other  frequency-dependent  input  information 
such  as  beam  pattern  parameters.  Beam  steering  angle  and 
array  orientation  are  also  entered  at  this  point.  The  character- 
istic function  program  reads  in  the  data  stored  in  the  geometry 
file  and  the  appropriate  transmission  loss  file.  The  program 
calculates  the  moan  and  variance  of  received  average  noise 
power  as  well  as  its  characteristic  function  and  stores  the 
results  on  a file  to  be  used  as  input  to  the  final  computer 
program . 
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The  final  program  computes  the  probability  density  function 
of  Y using  the  numerical  algorithm  described  in  Section  3 .2  of 
BBN  Report  3b'j3.  The  mean,  variance  and  probability  density 
function  for  Y are  printed  out  as  the  final  result. 

The  structure  of  the  program  provides  several  advantages. 
The  geometric  data  need  be  computed  only  once  for  a given 
sensor  location  and  route  structure.  Since  it  is  independent 
of  frequency  and  beam  steering  angle,  the  same  geometry  data 
can  be  used  many  times  by  the  characteristic  function  program. 
Furthermore,  the  transmission  loss  data  need  be  entered  only 
once  for  a particular  sensor  location. 

1 .2  Report  Overview 

The  remainder  of  this  report  documents  the  individual 
computer  programs  shown  in  Figure  1.  Section  J describes 
Program  GEOM  which  computes  the  geometry  data.  Also  described 
are  the  structure  of  the  input  data  file  of  sensor,  route  and 
ship  traffic  information;  and  the  structure  of  the  geometry 
data  file.  Section  3 describes  Program  CFUNK  which  computes 
the  characteristic  function  4>y(<*0  . The  structure  of  the  input 
and  output  files  needed  for  tills  program  are  also  described  in 
Section  3.  Section  4 describes  Program  DENS,  which  computes 
the  probability  density  function  Y and  prints  out  the  results. 

The  documentation  for  each  program  consists  of  a brief 
description  of  the  main  program,  the  structure  of  input  and 
output  files,  a glossary  of  variable  names,  flowcharts  and 
listings  of  the  program. 
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2.  PROGRAM  GEOM 

In  Section  2.2  of  BBN  Report  3653  were  derived  algebraic 
equations  to  be  solved  for  certain  geometric  parameters  needed 
in  the  computation  of  the  characteristic  function  of  received 
noise  'tv(u)).  These  geometric  parameters  depend  only  on  the 
position  of  the  sensor  and  the  position  of  shipping  routes  on  the 
earth's  surface.  They  can  be  broken  up  into  two  categories. 

In  the  first  category,  the  parameters  depend  on  the  variables  of 
integration  in  the  calculation  of  In  the  second  category 

they  do  not . This  memorandum  documents  a FORTRAN  computer 
program  called  Program  GROM  which  computes  the  geometric  para- 
meters of  the  second  category. 

Program  GROM  is  structured  to  compute  the  geometric 
parameters'  for  each  individual  segment  of  a shipping  route.  It 
is  assumed  that  a route  will  have  from  one  up  to  a maximum  of 
five  segments.  The  program  assumes  a maximum  of  ten  routes. 

The  program  reads  input  data  from  a previously  created  disk  file 
labelled  BATIN  and  stores  the  results  of  the  computations  on 
a disk  file  labelled  GEO.  The  results  are  also  printed  at  the 
computer  terminal  so  that  a hard  copy  is  available  for  documen- 
tation . 

2.1  Input  Data  File 

The  input  data  for  Program  GEOM  is  read  from  file  DATIN'. 
This  file  contains  the  identification  number  for  the  particular 
case  of  scnsoi  position  and  shipping  routes  to  be  processed;  the 
sensor's  position;  the  number  of  routes  and  the  appropriate  in- 
formation for  each  route  by  segments.  Table  1 gives  the  struc- 
ture of  the  file  DAT1N.  Each  line  of  the  table  corresponds  to 
one  logical  record  on  the  file. 

The  following  convention  was  assumed  for  the  position  of  the  sensor 
and  of  the  route  segment  end-points.  If  the  latitude  is  north, 
then  the  degrees,  minutes  and  seconds  that  specify  the  latitude 
are  all  positive.  If  the  latitude  is  south,  then  the  degrees, 
minutes  and  seconds  arc  entered  as  negative  numbers.  If  the 
longitude  is  west,  then  enter  the  degrees,  minutes  and  seconds 
as  positive  numbers.  If  the  longitude  is  east,  enter  them  all 
as  negative  numbers.  It  is  also  assumed  that  the  segments  are 
ordered  from  left  to  right  along  a route. 


2-1 


Report  3654 


Bolt  Beranek  and  Newman  Inc. 


WB 


TABLE  1 STRUCTURE  OF  INPUT  FILE  LATIN 

(Note:  Each  Line  Represent  One  Logical  Record) 


DESCRIPTION  FORMAT 

Identification  Number  (13) 


Sensor  Latitude,  Sensor  Longitude  (6F6.1) 

(in  degrees,  minutes,  seconds) 

Number  of  Routes  (13) 


Number  of  Segments  on  This  Route 

Latitude  and  longitude  of  lef t -end-point 
of  segment  (in  deg.,  min.,  see.) 

Width  of  route  segment  at  left 
end-point  (in  degrees) 

Latitude  and  longitude  of  right 
end-point  (in  deg.,  min.,  sec.) 

Width  of  route  at  right  end-point,  maximum 
width  of  segment  (both  in  degrees) 

Large,  small  and  fishing  ship  densities 
on  this  segment  (ships/n.  mile) 


etc . 


(13) 

(6F6 . l ) 

(F8.»l) 

(6f6 . 1) 

(2F8JI) 

(3F10.7) 
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2.2  Program  Description 

Figure  2 1.;  a flowchart  for  Program  GEOM.  The  program 
first  Initializes  constants  and  arrays;  defines  a function  KADS 
to  convert  angles  in  degrees,  minutes  and  seconds  to  radians; 
and  rewinds  the  input  and  output  files.  Next,  the  identification 
number  is  read  from  files  DATIN  and  also  the  sensor  position. 

The  sensor  latitude  and  longitude  are  convert'  d to  radio  The 

angle  between  the  north  pole  and  sensor  is  computed  as: 

_ - - latitude  of  sensor 

0 2 

The  angle  between  the  Greenwich  merridian  and  the  sensor  merridian 
is  computed  as: 

W =|sensor  longitude,  if  positive 
( 2 + longitude,  if  negative 

The  program  next  reads  from  file  DATIN  the  number  of 
routes  and  enters  a DO  LOOP  to  process  information  one  route 
at  a time.  The  first  statement  in  the  DO  LOOP  for  the  routes 
reads  in  the  number  of  segments  on  this  particular  route  and 
then  enters  a DO  LOOP  on  the  segments  of  this  route.  In  the 
DO  LOOP  on  the  segments,  the  program  reads  in  the  data  for  each 
segment  of  the  route  and  computes  the  geometric  parameters 
for  that  segment.  The  segments  are  ordered  from  left  to  right 
along  the  route.  The  geometric  parameters  and  shipping  density 
information  are  stored  in  a three  dimentional  array  G.  The  ele- 
ments of  array  G arc  defined  in  the  glossary  below.  The  com- 
putations performed  in  this  part  of  the  program  are  as  follows. 

First  the  angles  of  a north  polar  triangle  formed  by 
the  sensor,  route  end-point  and  north  pole  are  computed.  Figure 
3 shows  this  triangle  for  a general  end-point R.  . The  side  c.  is 
computed  from;  1 1 

c.  = - - latitude  of  end-point 

1 2 

The  angle  N^  is  given  by 


where 


longitude  of  end-point  if  positive 
2u  + longitude  of  end-point  if  negative 


2-3 


Report  3654 


Bolt  Beranek  and  Newman  Inc 


Figure  2 Flowchart  for  Program  GEOM 
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The  side  s^  of  the  polar  triangle  is  found  from 

cos  s.  = cos  c,  cos  c + sin  c,  sin  c cos  ft, 

1 1 o 1 o 1 

and  then  taking  the  arc  cosine.  Next  sin  s.  and  tan  - s.  are 

1 2 1 

computed.  The  angle  C.  of  the  polar  triangle  is  found  from  the 
arc  sine  of 


sin  = sin  sin  -f  sin  s^ 

The  azimuth  of  point  is  computed  as 


Z1  = 


Wi  1 


W 


= Pit  - C, 


Wi  > 


W 


The  triangles  formed  by  the  route  segment  end-points 
and  the  sensor  are  shown  jn  Figure  Jj . The  internal  angles  of  the 
route  segment  triangles  are  computed  as 


The  side  1^  of  the  route  segment  triangle  Is  found  from  taking 
tne  arc  cosine  of 


cos  1^  - cos 


s^  cot 


‘i+i + sln  5iti coi 


Then  the  length  of  the  route  segment  in  nautical  miles  is 
computed  from 


Li  h 

where  p is  the  radius  of  the  earth  in  nautical  miles. 

The  internal  ^ngle  F of  the  route  segment  triangle  Is 
computed  from 

sin  F.  = sin  s^  + ^ sin  sin  1^ 

and  then  taking  the  arc  sine. 
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V observation  point 


1' ! pure  4 Route  Segments  and  Triangles 
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The  last  computations  in  the  DO  LOOP  on  the  sop, merits 


is  to  compute  the  two  route  width  coefficients  b,  and  c.  for 
this  segment:  1 


1 

■ 

, "1 

h = 2 fc“ 

°i  * *i 

m^  - w^  + 

1 

*— 1 

4. 

•H 

£ 

1 

•H 

•H 

1 

2 r 

°i  ' *i 


2 r 

(mi  - vt^)  + (m1  - wi  + l^  + “ wi  ^ CmA  - wi  + 1) 


whei*e  w is  the  route  width  at  R^ 


rtu  is  the  maximum  width  of  the  segment  + 1 


After  the  computations  for  each  segment  of  a route 
have  been  completed,  the  program  returns  to  the  outer  DO  LOOT 
to  process  another  route.  Once  all  routes  have  been  processed, 
the  program  enters  the  section  of  code  which  writes  the  results 
to  the  output  file  GKO.  The  structure  of  the  output  file  GEO 
will  be  discussed  below. 


In  the  final  section  of  Program  GKOM,  the  results  of 
the  calculations  are  printed  at  the  terminal.  Figure  5 shows 
an  example  of  some  typical  printout  at  the  terminal . The 
first  thing  printed  is  a header  which  includes  the  identifica- 
tion number  and  the  number  of  routes.  Then  the  route  and  segment 
are  identified  and  the  thirteen  parameters  for  that  route  segment 
are  printed.  This  corresponds  to  one  row  of  the  array  G. 


2.3  Output  Data  File  GEO 


The  output  from  Program  GKOM  to  be  used  in  further 
computations  is  stored  on  disk  file  GEO.  File  GKO  is  one  of 
the  input  files  for  the  program  that  computes  the  characteristic 
function  of  received  noise,  Program  CFUNK . This  file  contains 
the  identification  number,  number  of  routes,  an  array  NS  containing 
the  number  of  segments  on  each  route,  and  the  array  G of  geometric 
parameters  and  ship  densities.  Table  2 gives  the  structure  of 
file  GKO. 
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2.4  Glossary  and  Program  Listing 


The  following  glossary  contains  definitions  of  the 
FORTRAN  variable  names  used  In  Program  GROM.  The  names  are 
presented  in  alphabetical  order. 

AN  = polar  angle  between  sensor  and  route  segment  end- 

point. measured  form  north  pole  (in  radians) 


C 


side  of  the  north  polar  triangle  formed  by  north 
pole  and  route  segment  end-point  (in  radians) 


CO 


side  of  north  polar  triangle  formed  by  north  pole 
and  sensor  location  (in  radians)  c 

o 


D1 

D2 

G ( I , J , K ) 


degrees  of  latitude  of  segment  end-point 
degrees  of  longitude  of  segment  end-point 

an  array  dimensioned  (5x13x10)  containing  the  computed 
geometric  parameters  and  ship  densities  for  each 

segment  of  each  route.  The  index  K indicates  the 
route  while  the  index  1 indicates  the  segment. 

For  a given  segment  on  a given  route  the  elements 
in  the  lih  row  are: 


G ( I , 1 , K) 
G(I,2,K) 


0(1, 3, K) 


G ( I , *1 , K ) 

G ( 1 , 5 , K ) 
G ( 1 , 6 , K ) 
G ( 1 , 7 , K ) 

G(I,8,K) 


route  segment  left  end-point  width  (radians) 

b^,  a parameter  used  to  model  the  route  segment 
width  envelope. 

e^, a parameter  used  in  the  model  of  the  route  seg- 
ment width  envelope. 

sin  Sp  sine  of  the  earth  centered  angle  s^  between 
the  sensor  and  a segment  left  end-point. 

cos  Sj, 

tan  - s. 

2 1 

F^  an  interior  angle  of  the  triangle  formed  by  the 
segment  and  the  sensor  (in  radians) 

the  azimuth  of  the  left  end-point  of  the  segment 
(in  radians) 


2-11 


Report  3654 


Bolt  Bcranek  and  Newman  Inc. 


G(I,9,K)  « 


Ij,  an  interior  angle  of  the  triangle  formed  by  the 
segment  and  the  sensor  (in  radians) 


0(1,10, K)  = length  of  the  segment  in  nautical  miles  L; 


0(1, 11, K)  » 


G ( 1 , 12 , K)  « 


0(1,13, K)  » 


IDNUM 


density  of  large  merchant  ships  on  this  segment 
(nhips/n.  mile) 

density  of  small  merchant  ships  on  this  segment 
(ships/n.  mile) 

density  of  fishing  vessels  on  this  segment  (ships/ 
n.  mile) 

an  ident i f ieai ton  number  for  the  particular  sensor/ 
route  geometry 

number  of  routes  (up  to  ten) 

an  array  dimensioned  10,  each  element  of  this  array 
contains  the  number  segments  on  a route. 


TEMP 
TEMP  2 
TEMP  3 


= radius  of  the  earth  in  nautical  miles 

= an  array  dimensioned  3 which  contains  the  sensor 
latitude  as: 

SLA  (1)  - degrees 
SEA  (2)  = minutes 
SLA  (3)  = seconds 

= an  array  dimensioned  3 which  contains  the  sensor 
latitude  as: 

SL0  (1)  = degrees 
SL0  (2)  = minutes 
SL0  (3)  3 seconds 

= seconds  of  latitude  of  segment  end-point 
« seconds  of  longitude  of  segment  end-point 


temporary  storage  locations 


longitude  of  end-point  in  radians 
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WIT  = width  of  route  segment  at  left  end-point 

WO  = longitude  of  sensor  in  radians 

XI  = minutes  of  latitude  of  segment  end-point 

^2  = minutes  of  longitude  of  segment  end-point 

7.  - temporary  storage  location 

The  listing  of  the  program  is  as  follows. 


PROGRAM  GFOM ( OUTPUT , PATTN, GFO, TA PE 2=DATI N , TA PE6 =GEO) 

C PROGRAM  TO  COMPUTE  THE  GEOMETRIC  COEFFICIENTS  NEEDED  IN 

C AMBIENT  NOTSE  MODEL.  WRITTEN  15  F FB . 78  BY  ZFSKINP. 

C DATIN  = FILE  OF  INPUT  SENSOR  AND  ROUTE  DATA. 

C CEO  = FILE  OF  GEOMETRIC  PARAMETERS  TO  PE  USED  AS 
C INPUT  TO  KORE  AND  LINE  COUNT  PROGRAMS. 

C 

DIMENSION  G(5, 13, 10),NS(10),SLA(3),SL0(3) 

DATA  G/650«0. 0/ , NS/1  0*0/ 

RADS(D,XM,S)=(1 .745329252E-02)*(D+(XM/60.)+ (S/3600. ) ) 
RHO=3H 37. 746771 
PI  = 3.  141592654 
REWIND  2 

REWIND  6 

C READ  IN  IDENTIFICATION  NUMBER 
R E A D ( 2 , 1 r' ) IDNUM 
C READ  SENSOR  POSITION 

READ(2, 10)(SLA(K) ,K=1, 3), (SLO(J) , J=1,3) 

10  FORM AT ( 6F  6.1) 

CO  = R ADS  (SLA  ( 1 ) , SLA  ( 2 ) , SLA  ( 3 ) ) 

CO=(PI/2. )-CO 

WO=RADS (SLO( 1 ) ,SLO(2) ,SL0(3) ) 

IF(WO. LT. 0. 0)  WO= (2. * P I )+WO 
T1=C0S(CC) 

T2=S IN ( CO) 

C READ  IN  NUMBER  OF  ROUTES 
READ(?.1S)  NR 
15  FORMATtI 3) 

C DO  LOOP  ON  ROUTES  *«**• 

DO  1010  K= 1 , NR 

C READ  IN  NUMBER  OF  SEGMENTS  ON  THIS  ROUTE 
READ(2 , 15)  N S (K ) 

NUMS  = NS(K) 

C DO  LOOP  ON  SEGMENTS  

DO  1000  J=1 , NUMS 
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DU  It  DUItmt-K  dliu  Numimn  lilt 


C READ  IN  SEGMENT  DATA  FOR  LEFT  END-POINT 
READ(2,  1 0)  D1,X1,S1,  D2,X?,S2 
C=(PI/?. )-RAPS(D1,X1,S1) 

W=RADS(D2,X2,S2) 

IF(W.LT.O.O)  W=(2.«PI)+W 
AN=ABS(W-WO) 

T 3 = C0S ( 0 ) 

T*=S1N  (C) 

G(J,5,  K)  = T^«T  WTU»T?«COS(AN) 

TEM  PrA  COS (G  ( J . G , K ) ) 

G ( J , « , K ) = S I N ( T EM  P ) 

IF(G ( J , 5 , K) . EO. - 1 .0)  GO  TO  5 

C(Jl6,K)=G(Jt^,K)/(1.+G(J,5,K)) 

GO  TO  6 

5 G( J , 6, K)= 1 . OE+99 

6 TEMP=TH*SIN(AN)/G(J , ^ , K ) 

G(J  , 8,  K)sASIN(TEMP) 

IF(G(J , 8 , K ) .LT.O.O)  G(J,8,K)=G(J,8,K)+(PI/2.) 
IF(W.GT.WO)  G(J,8,K)=(2.*PI)-G(J,8,K) 

C READ  IN  LEFT  END-POINT  ROUTE  WIDTH 
READ(2, 16)  WIT 
G(J, 1.K)s RAPS (WIT, 0.0, 0.0) 

C READ  IN  SEGMENT  RIGHT  END-POINT  POSITION 

READ(2, 10)D1,X1,S1,P?,X2tS2 
C=(PI/?.)-RADS(D1,X1,S1) 

W=RADS(D2,X2,S2) 

IF(W. LT.O.O)  W= (2. *P I )+W 
AN=ABS( W-WO) 

T3=C0S(C) 

TH=SIN(C) 

TEM P=T  3*T 1+T u«T?*COS ( AN) 

T EM  P2 = A COS  ( T EM  P ) 

TEMP3=S  IN  (T  F.M  P2) 

TEM  P?  =T ») »S  IN  ( A N ) /TEM  P3 
Z =A  SI  N (TEM  P2  ) 

IF ( Z . LT.O.O)  Z=Z+ (PI/2. ) 

IF(W.GT.WO)  Z= (2. * P I )-Z 
C COMPUTE  ANGLE  I 

G ( J , 9, K)  = ARS(Z-G(J , 8, K) ) 

C COMPUTE  LENGTH  OF  SEGMENT  L 

TEMP2=G ( J , 5,  K)«TEMP+G(J  ,4,K)  «TEMP?*COS(G(J  , 9,K)  ) 
TEM  P=A  COS  (TEMP?) 

G(J  , 10,  K)  = TEMP»RHO 
C COMPUTE  ANGLE  F 

TEMP?=T EMP3*SIN(G(J.9,K))/SIN(TEMP) 

G ( J , 7 , K ) = A S I N ( T EM  P? ) 

IF(G(J,7,K)  .LT.O.O)  G ( J , 7 , K ) = G ( J , 7 , K > + (P  T /.? . ) 

C READ  IN  RIGHT  END-POINT  WIDTH  AND  MAX  WIDTH  OF  SEGMENT 
R F. A D ( ? , 16)  W1D.WMAX 
16  FORM AT (2F 8. ^ ) 

WI D=R  APS ( W I D, 0.  ,0. ) 
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WMAX=R  ADS !WMAX, 0.  , P. ) 

C COMPUTE  ROUTE  WIDTH  COEFFICIENTS  FOR  THIS  SEGMENT 
T3=WMAX-G( J , 1 , K) 

T4=WMAX-WTD 

G(J  , 2,  K)=  !2./TEMP)  «(T?+S0RT!T?*T4) ) 

G(J  , 3»  K)  = (T,  + T4+SQRT(T'*T4>  )/(TFMP*TFMP) 

C READ  IN  SHIPPING  DENSITIES  FOR  SEGMENT 
READ(2,20)!G!J , L,K) ,L:1 1 , 13) 

20  FORMAT ( 3F 1 0. 7 ) 

1000  CONTINUE 
1010  CONTINUE 

C «•*»  OUTPUT  RESULTS  TO  FILE  GEO  «««* 

WRITE (6, 100)  IDNUM.NR 
100  FORMAT!! X, 13, ?X, 13) 

WRITE ( 6 , 1 10  ) NS 
110  FORMATC1 X, 1013) 

DO  115  K= 1 , NR 
NUMSrNS(K) 

DO  118  1=1 , NUMS 

WRITE(6, 120)(G(I , J,K) , J=1 , 13) 

120  FORMATd  X.4E20.  14) 

118  CONTINUE 
115  CONTINUE 

C PRINT  RESULTS  AT  THE  TERMINAL  

PRINT  200.IDNUM 

200  F0RMAT(1X///1X, 2 4H OUTPUT  FROM  PROGRAM  GFOM , 5X , 7HIDNUM 

PRINT  210, NR 

210  F0RMAT(1X/1X, 1 9H NUMBER  OF  ROUTES  = ,13///) 

DO  201 0 K = 1 , NR 
NUMS=NS(K ) 

DO  2000  Nr  1, NUMS 
PRINT  220,  K,  N 

220  FORMAT! 1X/1X, 6 H ROUTE  , 1 7,  ^X , 8HSEGMFNT  ,13) 

PRINT  230, G(N , 1 , K) ,G(N , 2, K) ,G(N , 3, K) 

230  F0RMAT(1X//1X,7HWTDTH=  , E 1 0.4 , 3X , 3HD=  fE10.4,3X, 

+ 3HE=  , E 1 0. 4 ) 

PRINT  240, (G (N , L, K) ,L=4,6) 

240  FORMAT! 1 X//1 X, 7HSIN  S=  , E 1 0. 4 , 3X , 7HCOS  S=  .E10.4.3X, 

♦ 10KTAN  1/2S=  , El 0.4) 

PRINT  250 , (G (N , M, K) , M = 7 , 10) 

250  FORMAT! 1 X//1 X, 3HF=  , E 10. 4 , 3X, 3HZ  = ,E10.4,3X, 

♦ 3HI=  ,E10.4, 3X, 3HL=  , E 1 0 . 4 ) 

PRINT  260,(G(N,J,K),J=11,  13) 

260  FORMAT! 1 X//1 X, 1 4HSH1 P DENS  IT IES/5X , 8HLA RGE  = ,E10.4, 

♦ 3X  , 8HSM ALL  = , E 10. 4 , 3X, 7HFISH  = ,E10.4)  - 
2000  CONTINUE 

2010  CONTINUE 
C 

REWIND  2 
REWIND  6 
STOP 
END 
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PROGRAM  CFUNK. 

The  algorithm  for  the  computation  of  the  characteristic 


function  of  beam  noise  power 
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to  implement  this  algorithm. 
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is  described  in  Section  2, 
propram  called  CHUNK  was  writ! 

Propram  CFUNK  computes  the 
part  of  the  characteristic 
function  of  beam  noise  power.  These  results,  alone,  with  other 
pertinent  information,  is  stored  in  a file,  called  PHIX, 
for  use*  by  Program  DENS . Program  DENS  computes  th  probability 
density  function.  It  is  described  in  Section  of  this  report. 

3.1  Input  Files 
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These  three  input  Piles  for  Program 
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3.1.2  Transmission  Loss  File 
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example,  transmission  loss  at  10  Herts  is  stored 
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3.1.3  File  SOURCE 

The  remaining  input  information  that  is  necessary  for 
program  CFUNK  is  read  in  from  a file  called  SOURCE.  This  file 
contains  frequency-dependent  input  information  and  beam  pattern 
parameters . 

Table  4 gives  the  structure  of  input  data  file  SOURCE. 

The  file  contains  the  mean  and  variance  of  the  total  radiated 
noise  for  large  merchant,  small  merchant  and  fishing  vessels 
at  the  frequency  of  interest.  This  is  stored  in  a 3 by  2 
dimensional  array  labeled  S.  The  center  frequency  and  frequency 
band  in  Hz  are  the  next  logical  record  on  the  file.  The  last 
logical  record  contains  the  beam  pattern  parameters.  They 
are  the  azimuth  of  array  broadside,  BAZ;  beam  steering  angle, 
BSTER;  a beam  pattern  constant  depending  on  frequency,  ALPHA; 
and  the  main  lobe  beam  width,  BWIDE.  All  angles  are  in  degrees. 
Figure  6 shows  the  definition  of  the  angles  associated  with 
the  array  pattern. 

3.2  Main  Program 

The  main  program  of  CFUNK  performs  the  following  functions: 
1)  initializing  arrays  and  set  program  constants;  2)  reads  in 
input  data;  3)  computes  initial  parameters  for  numerical  inte- 
gration; 4)  computes  the  probability  of  no  noise;  5)  computes 
the  mean  and  variance  of  the  averaged  noise  power  at  the  output 
of  the  beamformer;  6)  computes  the  real  part  of  the  character- 
istic function  of  averaged  noise  power;  and  7)  writes  the  results 
to  the  output  file  PHIX. 

The  following  subprograms  are  called  by  the  main  program: 

• Z(R,D)  - returns  a value  of  transmission  loss 
for  a given  range  and  angle . 

• BETA  (Q,  QO)  - returns  a value  of  the  probability 
density  function  of  transverse  ship  position 
across  a route  segment. 


• BEAM  - a subroutine  to  compute  the  normalized 
beam  pattern  for  a given  angle . 

• GAMA  - a subroutine  to  compute  the  real  and 
imaginary  parts  of  the  characteristic  function 
of  the  gamma  probability  density  function. 

• INCLU  - a utility  function  subprogram  which 
returns  a value  of  one  if  an  angle  C is  included 
between  two  angles  A and  B. 


These  subprograms  are  described  in  Section  3.5. 
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3.2.1  Common  Storage 

There  Is  one  block  of  common  storage  area  labeled  XLOSS. 
XLOSS  contains  the  array  TL  of  transmission  loss  data  and  is 
common  to  the  main  program  and  the  function  subprogram  Z. 

3.2.2  Flowchart  and  Description 

Figure  7 presents  the  flowchart,  of  the  main  program.  At  the 
start  of  the  program,  the  arrays  (i , NS  and  PHI  are  Initialised  to 
zero.  Program  constants  are  set  as  follows: 


HI  10 

3^37. 746771 

(radius  of 

the  earth  in 

nautical 

m i les ) 

PI 

3.1411192654 . 

TPI 

2 n 

DELR  = 

0.0029088821 

radians . 

(range  increment 

step  size 

equivalent  to  10 

nautical 

mi les ) 

DIO 

0.1745329252 

radians . 

(equivalent  to  10°) 

The  following  statement  functions  are  defined  at  the  start 


the 

program : 

• 

RADS  - 

converts  angles  in  degrees,  minutes  and 
seconds  to  radians. 

• 

WIDTH  - 

computes  an  approximation  to  route  segment 
width 

• 

XJ 

computes  the  value  of  the  Jacobian  of  the 
transformation  between  sensor  coordinates 
and  route  coordinates. 

The  program  next  reads  in  the  file  GEO  containing  the 
geometric  parameters  computed  in  Program  GliOM  and  stores 
them  in  array  G.  Next,  the  appropriate  transmission  loss 
data  is  read  from  file  and  stored  in  array  TL.  The  first 
row  of  TL  is  converted  from  degrees  to  radians.  Finally, 
the  file  SOURCE  is  read.  This  completes  the  Input  section 
of  the  program. 
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IN  IT IAI  I 71 
511  PROGRAM 


M ARRAYS  AND 
KAN  CONSTANTS. 

I NO  mis 


riad  muwmtt)  raomiRv  data 
I KOM  | at  GlOMl  IRY  . 

STORl  IN  ARRAY  G 

I 


R(AO  IN  TRANSMISSION  LOSS 
DATA  f ROM  I III . 

STORl  IN  ARRAY  Tl  j 


I C ON VI RT  I 
I I ROM  DIC.H 


RT  MHST  ROW  Of  Tl 
DIGRUS  TO  RADIANS. 


RIAD  IN  DA  I A ON 
Till  SOUKCt.  i 


C Al  ClM  All  I'ARAMIURS  01  GAMMA  PROOApItUY  DINSITY 
f UNCI  ION  IRON  MIAN  ANO  VARIANCE  01  SOURCTS. 

CONVIRT  ANGUS  Of  M NSOR  IRON  Of  GRl  l S TO  RADIANS. 


INITIAI  CALCULATIONS 

COMPUIT  PACMUAM  AND  T HI  TOUR  ANGl  I RIG  I ONS  01  INTIGRAT  ION. 
COMT’UT  l STII1  SI/IS  I OR  ANGl  I INTIGRAT  ION  ANO  NUMBLR  OF 
ANO 1 1 IT  I NATIONS. 

T T ST  FOR  TNO  MKI  CONDITION  SIT  I'ARAMIURS  I OR  ANGl  l LOOP. 


COMPUU  1 Ml  PROPAD  1 1 1 T Y 
01  NO  NO  I SI  . Dili. 


tSIT  XMIAN  - 0 1 

XVAR  - 0 


DO  LOOP  ON  ANGl l 
I A - T , NT 


SCltCT  WHICH  ANGl  l INTIGRAT  ION  Sill’  SlZl  TO  USt  PI  PL  ND I NG 
VmiHIR  ANGl  ( IS  IN  MAIN  HI  AM,  PAC  NR  l AM  OR  SIDUOPIS. 

CONVIRT  ANGl  I INTO  CQIVALCN1  Dl  A/l  PAT  URN  ANGL  l I ROM 
PROADSIDl.  UT  no  . o.o 

Til  - 0.0 


f CALL  Pf AM 

i hi  turns  hi  am  pat  urn  vaiui 

\ I OR  GIVI N ANGl l / 


00  100P  ON  ROUTT S 
K - \ . NR 


00  100P  ON  ROUTT  SfC.MINlS 


/ IS 

THIS  ANGl  IN.  NO 
UNCI  01*  l 0 IN  THIS?— 


FIGURE 


Flowchart  of  Program  CFUNC. 
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inc  jju  l t jyjt 


du  i i Deru lie k dim  ncwumii  me 


The  program  next  enters  a section  of  code  which  performs 
the  Initial  calculations.  It  calculates  the  two  parameters 
for  a gamma  probability  density  function  from  the  means  and 
variances  of  the  sources  that  were  input  from  file  POUKdl  for 
large  merchant  ships,  small  merchant  ships  and  fishing  vessels. 
That  is,  each  of  the  three  sources  Is  modeled  as  a random 
variable  that  is  gamma  distributed.  The  two  parameters  of 
the  gamma  probability  density  function  are  computed  as: 


and 


M. 


b+1  = -- 


V 


,th 


where  Mi  is  the  mean  of  the  j*"1  source  and 
The  program  then  converts  the  angles  assoc 
pattern  of  the  array  from  degrees  to  radians 


d 

i;it  bd 


is  its  variance 
with  the  beam 


In  the  final  part  of  the  initial  calculations,  t he  program 
computes  the  parameters  needed  later  in  the  program  for  tin-  DO 
loop  on  azimuthal  angle  from  the  angles  associated  with  the 
sensors  beam  pattern.  From  the  steering  direction,  bear. width 
and  array  broadside  angles,  the  program  finds  the  backbeam.  ;• 
also  tests  for  an  end  fire  condition.  If  an  end  fire  condition 
exists,  it  divides  the  angles  into  two  regions — main  beam 
and  side  lobe.  If  an  end  fire  condition  does  not  exist,  it 
divides  the  angles  into  four  regions--main  beam,  back  beam. 


ami  two  side  lobe  regions.  The  program  sets  the  angle  inte- 
gration step  size  equal  to  1/16  the  boamwidt h in  the  main  lobe 
and  back  lobe  regions.  In  the  side  lobe  regions  the  angle 
integration  step  size  is  set  equal  to  approximately  10  degrees. 
The  program  computes  the  total  number  of  angle  integrations  NT 
to  cover  the  complete  3t>0  degrees,  which  is  used  in  the  PO 
loop  on  angle.  This  is  the  end  of  the  initial  section  of  the 
program . 
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The  program  next  enters  a section  ol'  code  which  computes 
the  probability  of  no  noise,  with  variable  name  PEl.Z.  This 
term  was  derived  in  Section  2.1  of  BBN  Report  36r,>3,  and  is 
computed  as 


DELZ  = 


m n 

k.  . 
ij 

1=1  j = l 


i, 

i 


where 


is  the 
miles  and 
of  type  j 


length  of  the  i*’*1  route  in  nautical 
ki j is  the  average  number  of  ships 
per  nautical  mile  of  route  i. 


In  the  next  section  of  code,  the  mean  and  variance  of  the 
averaged  noise  power  at  the  output  of  a beamformer  is  computed. 
The  equations  were  derived  in  Section  3.1  of  BBN  Report  36^3. 

The  equations  are: 


m n 

n'x  = p T.  X k..Ms  f dD  / dr  | J ( r , D ) | f (q ) s (r  ,D ) . Bi'(D) 
i = l .1=1  J j R.  wi 


o 


2 

X 


m n 

p X X k,  . (m5  + o2  )f  dP  / dr  | J (r  ,D)  I fn  (q ) s2  ( r , B)  . BP 

i=l  .1  = 1 Sj  R.  Qi 


where 


m^,  is  the  average  of  the  mean-squared  pressure 
“ j of  sources  of  type  J 

o2  is  the  variance  of  the  mean-squared  pressure 
“"j  of  sources  of  type  ,|  . 


These  are  program  inputs  from  file  SOURCE. 


: IP) 
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i S 


lit  t ' 


Tit ' !1  * 


The  inner mo 

the  DO  loop  on  the  route  segments  to  a loop  on  the  range.  Thus 
the  program  evaluates  the  double  integrals  for  the  mean  anti 
variance  by  fix  inf.  the  angle  and  evaluating  the  integral  over 

in  v he 
, until 


range 

across 

t he 

route  segment  w .1  d * 

h t o 

compi 

l 1 0 

one 

t o rm 

angle 

int ogr 

a 1 . 

It  then  goes  on  to 

the 

next 

any 

tie, 

etc  . 

it  ha 
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et  ed 

the  calculation. 

After  c 

alculating  the  mean  and 
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noise 

power , 

t he 

program  re-init ial 
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cert ; 

i j n 

par 

amei  t 

that 

cent  rol 

t he 

numerical  int egrat 

ion  c 

md  ei 

it  er 

-s  l 

he  se 

of  code  that  computes  the  real  part  of  the  characteristic 
function  of  the  average  noise  power  v>').  The  real  part  c 
$ (m)  is  used  In  Program  OHMS,  described  in  Section  4 below 
to  compute  t lie  probability  density  function  of  average 
power . 


r ! 


.01 


\ 


The  argument  of  4>y  labeled  F in  the  program,  is  initialise 
to  aero.  The  number  of  points  at  which  the  real  part  of  4-y 
is  computed  is  set  to  NF  = 101,  with  a step  sice  of  DELS  = 0 
If  more  points  or  a different  spacing  foi- 
ls desired,  the  vnltu  of  ME  and  DKi.F  must 
program.  The  real  part  of  4 y is  computed 
In  the  first  step,  the  real  and  imaginary 
character!  st  ic  fund  ion  H'tw)  are  computed 
(See  Section  2 of  DBN  Report  3o‘>3' 


i. 


t he  argument  i 
be  changed  in 
in  two  basic  steps, 
parts  of  the  second 
from  the  equal  ion 


m r r j 

y(u-)  = p Z I I d | f (q)  2i  k 4* 
i=l  ' J Wi  1=1  ''  b- 


(u-  BP  2.)  dr  dp 


D r 


i 


where  4*  is  the  gamma  probability  density  function  model  for 

^ t 

*'  the  source  of  type  j . 
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This  computation  is  accomplished  by  using  rectangular 
integration  to  numerically  evaluate  the  double  inle.-al. 
This  is  done  in  a set  of  five  nested  loops.  The  outermost 
loop  being  on  the  angle  D,  while  the  innerr.  st  is  , the 
argument  of  the  characteristic  function  F. 


Once  the  second  characteristic  function  has  been  computed, 
the  program  computes  the  real  part  ol'  the  character istie  function 
from: 


Heal  ,(w)  = DEI, 


EXP  ^Rea  1 4'  ( w • COS  ^ 


Imaginary  41  (v 


a) 


for  the 
tat  ions 
the  real 
PHI  . 


101  values  of  the  argument.  This  concludes  the  compu- 
tes: Lons  of  Progri  n CFUNK . The  result ing  values  of 
pari  of  and  the  value  of  w are  stored  in  array 


In  tiie  terminal  section  of  the  program,  the  results  of 
Program  CHUNK  are  written  to  disk  file  PHIX,  which  is  described 
in  the  next  section.  The  input  and  output  files  are  then 
rewound,  and  the  programs  terminates  execution. 


3.3  Output  File 


The  output  from  Program  CFUNK  to  be  used  in  further  compu- 
tat  L< ns  s stored  on  disk  file  PHIX.  File  PHIX  is  the  Input 
file  for  Program  DENS  described  in  Section  4 below.  This 
file  contains  the  header  "RESULTS  FROM  PROGRAM  CFUNK"  as  the 
first  logical  record.  The  remainder  of  the  file  contains  the 
values  of  the  parameters  for  the  array;  the  total  number  of 


ships;  the  probability  of  no  noise; 
of  the  beam  noise  power;  the  number 


been  computed  for  the 
the  real  part  of  * 
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Y 


real 
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t he  mean  and  variance 
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5 gives  the 
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3.4  Glossary  and  Program  Listing 


The  following  glossary  contains  the  definitions  of  the 

FORTRAN  variable  names  used  in  Program  CFUNK.  The  names 

are  presented  in  alphabetical  order. 

ALPHA  - beam  pattern  constant  which  depends  on 
frequency.  Program  input. 

ANG  - azimuthal  angle  from  north  in  radians; 

initialized  to  EDGE  1 and  used  as  angle 
of  integration. 

BACK  - azimuth  of  back  beam. 

BAZ  - azimuth  of  array  broadside;  program  input. 

BP  - value  of  normalized  beam  pattern  returned 

from  Subroutine  BEAM. 

BSTER  - beam  steering  angle  from  north,  positive 

clockwise;  program  input. 

BWIDE  - beam  width  of  main  lobe  of  the  array; 
program  input . 

B1F  - (b  + 1)  parameter  for  gamma  probability  density 

function  model  of  fishing  vessel  radiated  noise, 
equal  to  mean  squared  divided  by  variance.  Used 
as  input  to  Subroutine  GAMMA. 

B1L  - (b  + 1)  parameter  for  gamma  probability  density 

function  model  of  large  merchant  ship  radiated 
noise . 

BIS  - (b  + 1)  parameter  for  gamma  probability  density 

function  model  of  small  merchant  ship  radiated 
noise . 

CP’  - c parameter  for  gamma  probability  density 

function  model  of  fishing  vessel  radiated 
noise,  equal  to  mean  divided  by  the  variance. 

CIF  - imaginary  part  of  the  characteristic  function 

for  fishing  vessel  radiated  noise  returned  by 
Subroutine  GAMA. 
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OIL 

CIS 

CL 

CS 

CKK 

CRL 

CRS 

D 

DANG 

DANGE 

DANG  1 

DANG  2 

DANG  3 


imaginary  part  of  the  characteristic  function 
for  large  merchant  ship  radiated  noise  returned 
by  Subroutine  GAMA. 

imaginary  part  of'  the  characteristic  function 
for  small  merchant  ship  radiated  noise  returned 
by  Subroutine  GAMA. 

c parameter  for  gamma  probability  density 
function  of  large  merchant  ship  radiated 
noise  . 

c parameter  for  gamma  probability  density 
function  of  small  merchant  ship  radiated 
noise . 

real  part  of  the  characteristic  function 
for  fishing  vessel  radiated  noise  returned 
by  Subroutine  GAMA. 

real  part  of  the  characteristic  function 
for  large  merchant  ship  radiated  noise 
returned  by  Subroutine  GAMA. 

real  part  of  the  characteristic  function 

for*  small  merchant  ship  radiated  noise  returned 

by  Subroutine  GAMA. 

deviation  angle  from  left-end  of  a segment, 
in  radians;  used  in  numerical  integration. 

step  size  in  deviation  angle  (AD)  to  be  used 
In  numerical  integration. 

deviation  angle  step  size  to  be  used  in  end- 
fire  case  . 

deviation  angle  step  size  in  the  main  lobe  of 
the  beam  pattern  ; set  to  l/lt<  of  the  beam 

width . 

deviation  angle  step  size  to  be  used  In  region 
between  main  lobe  and  back  lobe  going  clockwise. 

deviation  angle  step  size  to  be  used  in  the 
back  beam;  set  to  1/16  of  beam  width. 
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DANA  4 

DELF 

DELFQ 


DELE 

DELZ 
DIO 
EDGE  1 

EDGE  2 


deviation  angle  step  size  to  be  used  in 
region  between  back  lobe  and  main  lobe, 
going  clockwise. 

step  size  of  the  argument  of  characteristic 
function  (Aw)  of  beam  noise  power;  fixed  in 
program  to  0.01. 

frequency  band  of  received  noise  in  Hertz. 

integration  step  size  for  range  r set  in 
radians  equivalent  to  10  n.  miles  on  earth's 
surface . 

computed  probability  of  no  noise. 

a program  constant  in  radians  equal  to  10  degrees. 

azimuth  of  left-hand  edge  of  main  lobe  of 
beam  pattern  in  radians. 

azimuth  of  right-hand  edge  of  main  lobe  of 
beam  pattern  in  radians. 


EDGE  3 

EDGE  4 

F 

FREQ 
G ( I , J , K ) 

G ( I , 1 , K) 


azimuth  of  left-hand  edge  of  back  lobe  of 
beam  pattern  in  radians. 

azimuth  of  right-hand  edge  of  back  lobe  of 
beam  pattern  in  radians. 

argument  w of  the  characteristic  function  of 
beam  noise  power  <t»  (u>). 

center  frequency  of  band  of  received  noise 
in  Hertz. 

an  array  dimensioned  (5x13x10)  containing  the 
computed  geometric  parameters  and  ship  densities 
for  each  segment  of  each  route.  The  index  K 
indicates  the  route  while  the  index  I indicates 
the  segment.  For  a given  segment  on  a given 
route  the  elements  in  the  :iUi^  row  are: 

route  segment  left  end-point  width  (radians)  UL  . 
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0(1,2, K) 

b.  a parameter  used  to  model  the  route  segment 
width  envelope. 

0 (1 , 3 , K ) 

- 

e^  a parameter  used  in  the  model  of  the  route 
segment  width  envelope. 

G ( 1 , 4 , K ) 

- 

sin  Sj,  sine  of  the  earth  centered  angle  a. 
between  the  seimor  and  a segment  left  end-point  . 

G ( I , 5 , K ) 

- 

cos  s . 

0(1, 6, K) 

- 

tan  1/2  s . . 

0(1, 7, K) 

- 

an  interior  angle  of  the  triangle  formed  by 
the  segment  and  the  sensor  (in  radians). 

0 ( 1 , 8 , K ) 

- 

Z , the  azimuth  of  the  left  end-point,  of  the 
segment  (in  radians). 

0(1, 9, K) 

- 

1^  an  interior  angle  of  the  triangle  formed  by 
the  segment,  and  the  sensor  (In  radians). 

G ( I , 10 , K ) 

- 

length  of  the  segment  in  nautical  miles  I.,. 

0(1, 11, K) 

- 

density  of  .Large  merchant  ships  on  this  segment 
(shlps/n.  mile). 

G (I , 12, K) 

- 

density  of  small  merchant  ships  on  this  segment 
(ships/n.  mile). 

G(I,13,K) 

- 

density  of  fishing  vessels  on  this  segment  (ship 
n . mile). 

01 

- 

along, -route  position  variable  g measured 
from  left  end  of  a segment  . 

GIO 

- 

initial  value  of  01. 

1DNUM 

- 

an  identification  number  for  the  particular 
sensor/route  geometry . 

IKNDF 

- 

flag  for  end-fire  condition  of  beampattern 
where 

IENDF  = 0 Not  endfire  condition 
I1SNDF  = 1 Endfire  condition. 
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.1 RKLAU 


1 TEST 


IW 


N b' 


NR 

nrl:no 


flap  used  to  determine  how  and  when  to  reset 
ranp,c  loop  in  ranpe  intepration  where 

IRFLAG  = 0 then  1/2  width  of  the  route  is 
less  than  10  n.m. 

IRFLAG  = +1  In  upper  half  of  width  of  route 
sepment 

IRFLAG  = -1  In  lower  half  of  width  of  route 
segment . 

flap  to  determine  if  angle  AN(i  lies  in  a route 
sepment,  returned  from  function  INCLU,  where 

1TEST  = 0 Not  in  segment 
I TEST  = +1  In  sepment. 

index  of  DO  loop  on  argument  of  characteristic 
function  of  beam  noise  power. 

number  of  sample  points  at  which  characteristic 
function  of  beam  noise  power  is  computed;  set  to 
101  in  program. 

number  of  routes  (up  to  ten). 

number  of  angle  Iterations  for  numerical  inte- 
gration in  endfire  case. 


NR1 

number 

of 

angl  e 

NR2 

number 

of 

angle 

back  lobe 

iterations  in  main  lobe. 
Iterations  between  main  lobe  and 


NR3 

number 

NR 'l 

number 

and  inn 

■ of  angle  Iterations 

• of  angle  iterations 
in  lobe. 


in  back  lobe, 
between  back  lobe 


NS 


NT 

Nil  MS 


an  array  dimensioned  1.0,  each  element  of  this 
array  contains  the  number  segments  on  a route 

total  number  of  angle  iterations  in  the  angle 
DO  loop. 

temporary  storage  for  number  of  segments  of  a 
route . 
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I'll!  (.’,101  ) 


P 1 


w' 


K 


RAUF 


KHO 

HO 


O 

k' 


(3 


) 


SUM  (2  .101) 


TL( 301 , b ) 


Tl'  l 


TOtl  IF 


an  array  dimensioned  (2  x 301)  which  contain:-  tin 
computed  values  of  the  real  part  of  the  charactei 
istlc  function  i'^.(^)  and  w where: 

PHI(1,K)  = real  part  of  1,.  ( w ) 
rill  (2 , K)  = w 

This  array  is  output  to  file  PH1X. 

ii 


across  route  position  variable  q • initially 
ret  to  zero. 

range  from  sensor;  earth  centered  angle  in 
rad i ans . 


statement  function  to  convert  angles  in  degrees, 
minutes  and  seconds  to  radians. 

radius  of  the  earth  in  nautical  miles. 


range  from  sensor  to  the  center  of  a route 
segment  (earth  centered  angle  in  radians). 

an  array  dimensioned  (3  x 2)  containing  the 
values  of  the  mean  and  variance  of  total 
radiated  noise  for  the  three  sources,  where: 


S(l,l) 
S( 1 , 2 ) 
S(2,l) 
S ( 2 , 2 ) 

5(3,1) 

5(3,2) 


mean  of  large  merchang  ships, 
variance  of  large  merchant  ships, 
mean  of  small  merchant  ships, 
variance  of  small  merchant  ships, 
mean  of  fishing  vessels 
variance  of  fishing  vessels. 


an  array  dimensioned  (2  x 101)  used  to  contain 
the  running  sum  in  the  number! cal  integration 
over  range  part  of  the  calculations  of  the 
second  characteristic  fund  ion  where 


SUM(1,K)  = real  part 
SUM(2,K)  - imaginary  part. 

an  array  dimensioned  (3‘>1  x 5)  containing  the 
transmission  loss  table  read  into  the  program. 

2 it . 

Total  number  of  ships 
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TEMP  1 v 
TEMP  2 
TEMP  3 
T1  / 


T2 


T3 


T10 


Til 


Temporary  storage  locations 


V 


WI 

WIDTH 


X 


XJ 


XMEAN 


argument  for  source  characteristic  function 
used  in  calling  Subroutine  GAMA;  equal  to 
w • BP  • Z. . 

l 

one  half  the  width  of  a route  segment . 

statement  function  to  compute  an  approximation 
to  the  I’oute  width  between  route  segment  end 
points . 

angle  ANG  converted  to  beam  pattern  coordinates 
used  for  input  to  Subroutine  BEAM. 

statement  function  to  compute  the  absolute  value 
of  the  Jacobian  of  the  transformation,  where 

J = sin  in /cos  q .. 

computed  value  of  the  mean  of  received  beam 
noise  power. 


XS 


beam  steering  angle  with  respect  to  array 
broadside,  used  as  input  to  Subroutine  BEAM. 


XVAR 


computed  value  of  the  variance  of  received 
beam  noise  power. 


o o o o 
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The  follow  1 ng  is  a listing  of  the  FORTRAN  code  for 
the  main  program. 


PROGRAM  fl'l'KK  (OUT  TUT  , GKO,  TI.1PI!  , SOURCE , I’lll  )' , T A IT  IsGl'O, 
4 T A PE  2=T1. 1 0!) , TAPF  -USOURC  E,  TAPEGrPII  1 X) 

C 

c 

C PROGRAM  TO  COMPUTE  Til!'  C HA RACTE R 1ST  TO  FUNCTION  OF  RFCI'lVFp 


C NOISE  AT  THE  SENSOR.  WRITTEN  JUNE  1078  BY  K.  7ISK1ND. 

C 

c INPUT  KITES  

C CEO  = CONTAINS  ROUTE/SENSOR  CEOF'ETR  1C  PARAMETERS  A NO 
C SHIPPING  DENSITIES  COMPUTED  IN  GEOMETRY  PROGRAM 

C Tl.lOUr  CONTAINS  TRANSMISSION  TOSS  DATA 
C SOURCEr  CONTAINS  SOURCE  RADIATED  NOISE  DATA  , P I A M 
C PATTERN  PARAMETERS  A HP  OTHER  FREQUENCY  DEPENPENT 

C PARAMETERS. 

0 OUTPUT  FITE 

C PHIXr  CONTAINS  REAI.  PART  OF  CHARACTERISTIC  FUNCTION  AS 
C COMPUTE P BY  THIS  PROGRAM  AND  THE  MEAN  AND  VARIANCE. 

C USED  AS  INPUT  TO  THE  PR OBAB1 1.1  T Y DENSITY  PROGRAM. 


C 

C 

DIMENSION  G (S,  11,10),  NS(  10),  PH  UP,  1 Cl  ),*(?,?)  , SUM(P,  101 
COMMON /XTOSS/TT(  1 r>  1 , 5 ) 

DATA  0/650*0.  (V,  NS/  10*0/,  PH  I /P0?*0.  0/ 

C 

C STATEMENT  FUNCTIONS  AND  PROGRAM  CONSTANTS 

R A DS  ( D , XM  , S)  = 1 . 7 '1 5 IPOS  DE-OP*  (D«  (XM/60.  )<  (S/3600.  ) ) 

WIDTH  (XG  , P,  E,  XW)r  XK-.  (h*XG)+  (E*XG*XG) 

XJ  (R.1 , 01  ) --APS  (SIN  (R  I ) /COS  (01  ) ) 

DETRrP. 90888P1E-0 3 
R HO  = 3 A X 7 . 7A677  1 
PI =1.  1 A 1 5 9 P 6 5 A 
TPI=2.*PT 
D 1 0 = 0 . 17A5.3POP5? 


REWIND  FITES 
REWIND1 
REWINDS 
REWIND 3 
REWIND 6 
C 
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C READ  TH  TATA  FROM  FILE  CFO 
REAP(1,10)  1PNUM.NR 
10  FORMAT (IX, IT, OX, IT) 

K FA r> ( 1,1'))  NS 
15  FORMAK  IX,  1013) 

DO  ?0  K- 1 , NR 

NUMSrNS (K ) 

no  I - 1 , Kill's 

RKADC 1 , 30  > (C ( I , J,K),J  = 1,  13) 

TO  FORMAK  1 X.  REPO.  1 R ) 

?5  CONTINUE 

?0  CONTI NUF 

C 

C READ  IN  TRANM1SSI0N  LOSS  PATA 
READ(P,  *10)  T 1,  TO,  T T 

RO  FORMAK  1 X,  A 10,  A 10,  A5) 

R FAP  (?,50)((TL(K  1 , K?)  , KP=  1 , 5 ) , K 1 = 1 , T51  ) 

50  FORK  AT(  1 X , 5F.  1 !>.  8 ) 

C 

C CONVERT  FIRST  ROW  OF  TL  FROM  DEGREES  TO  RADIAL'S 
IX)  60  .1  = 1,5 

T L ( 1 , .1 ) = R A PS  ( T L ( 1 , J ) , 0 . ,0.) 

60  CONTINUE 

C 

C READ  IN  SOURCE  PATA  FILE 


IX)  70  K = 1 , 3 

R FA  D(  3,80)  S(K, 

1 ) ,S(K , ?) 

BO 

FORMAK  IX,  PESO. 

IT) 

70 

CONTINUE 
REAPH.BO)  FREQ 

, PEI.F Q 

REA  D( T , ps  ) P A 7. , 

PSTER, ALPHA , PW1PF 

B5 

FORM  AT (RE  IB.  10) 

C 

C COMPUTE  C ANP  IKI  GAMA  DENSITY  FUNCTION  PARAMETERS 
C FROM  SOURCE  MEANS  ANP  VARIANCES  FOR  LARGE,  SMALL  ANP 
C FISHING  VESSELS 

CL=S ( 1 , 1)/S(1 , ?) 

B11.  =CL*S  ( 1 , 1 ) 

CSrS(P,  1 )/S(P,P) 

U1G:CS*S(p,  1 ) 

CFrSn,  1 )/S( T, ?) 

B1F  = CF*S(3,  1 ) 

C 

C CONVERT  UA7  , PSTER  ANP  PW1PE  TO  RADIANS 
PA'/=R  APS  ( UA7 , 0.  ,0.  ) 

PST E h = R A PS  ( P ST E R , 0 . , 0 . ) 

PW  1 PE  = R A l\S  ( PW  ] PE  , 0 . ,0.  ) 

C 

c iMtmHKHMHiHinH*  END  OF  INPUT  SECTION  «*n<*»«**» 


A 


• • 
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INITIAL  CALCULATIONS 

COMPUTE  ANGLE  INFORMATION 
COMPUTE  AZIMUTH  OF  MAIN  LOBE  EDGES  El  AND  E2 

EDGE  1 rBSTER- (0. 5*BWIDE ) 

IF  ( E DGE  1 . LT . 0 . ) E DGE  1 =T PI  + E DGE  1 
F.DGE2=B  STER+0 . 5 * PW  I DE 
IF  ( E DGE  2.  GE.  TPI ) K DGE  ?-?.  DGE  2-T  PI 
C FIL'D  THE  BEAM  STEERING  ANGLE  WITH  RESPECT  TO  BROADSIDE 
C AND  THE  STEERING  ANGLE  FOR  SUBROUTINE  BEAM: 

TEMPI  =BSTER -B AZ 

I F ( R AZ  . GT  . BST E R ) T EM  P 1 =T  EM  P 1 +T  P I 
XSrTEKPI 

I F ( T EM  P 1 . GT  . PI  )XS  = XS-TPI 
C TEST  FOR  ENDFJRE  CONDITION 
I F.KDF =0 

I F ( X S . F.O.  PI.  OR.  XS.EG.-PI)IENDF=1 
C COMPUTE  BACK  BEAM  STEERING  ANGLE  AND  EDGES  E3  AND  EA 
BACKrPI-TEMPI 

IF(TEKP1.  GT.  PI  ) BAC K =B ACK+T PI 

EDCE3=BACK-0.5*BWTDE 

IF  (EDGE  '3.  LT. 0 )EDGE  3=E  DGE  3+TPI 

EDGE4=BACK+0.5*BWIDE 

IF  ( E DGE 11 . GE . TPI )E  DGE  A = E DGE4-TPI 

C INITIALIZE  ANGLE  TO  EDGE  1 AND  SET  ANGLE  REGION  1 PARAMETERS 
A NG  = E DGE 1 
DANG 1=BWIDE/1 6.0 
N R 1 = 1 6 

C COMPUTE  PARAMETERS  FORR  THE  REMAINING  REGIONS 
IF ( I ENDF  . EC.  1 )GOTO  115 
C REGION  TWO 

T E M P 2 = E DG  E 3 - E DG  E 2 

IF  (TEM  P2.  LT.  0. 0 )TEM  P2=TEM  P2+TPI 

NR2=I NT (TEM  P2/D 10) 

DA  NG2=T EM P2/FLOAT (NR2) 

C REGION  THREE 
NR3=N  R 1 
DANG3=DANG1 
C REGION  FOUR 

TEMP3  =E  DGE 1 -E  DGEA 

IF(TEMP3.  LT.  0.  0)TEMF3=TEMP3+TPI 

NRA=INT(TEMP3/D10) 

DA  NG  U =T  EM  P 3 /FLOAT  ( N R 'I ) 

C COMPUTE  PARAMETERS  FOR  A NG  DO  LOOP 
NT  = N R 1+NR2+NR  3+NR A 
GOTO  118 
C ENDF 1 RE  CASE 
115  T EM  P 3 =T  P I -B  W I DE 

NRENDr I NT (TEM  P3/D 1 0 ) 

DANGE=TEMP3/FLOAT(NREND) 

NT  = NR  1+NREND 
118  CONTINUE 

£*#«***#*##»*«##*# end  Q]r  INITIAL  SECTION*4*  ********** * 

• • 
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C COMPUTE  THJ-  PR  OB ABILITY  OF  NO  NOISE , TERM  DEL? 

TSUI  P=0 . 0 
DO  IPO  K = 1 , NR 
NUMS=NS(K  ) 

PO  12B  1 = 1,  HUMS 

TSII1  PrTSMI  P+G  ( I , 1 0,  K ) * (G  ( I , 1 1 , K ) *G( T , I ?, K) 4 G ( 1 , 1 3,  K ) ) 
IPr>  CONTINUE 
IPO  CONTINUE 

DELZ  = EXP( -TSHI P) 


SECTION  TO  COMPUTE  MEAN  AMP  VARIANCE 
XMFAH=0. 0 
XVAR  = 0. 0 
PO  POO  I A = 1 , NT 
T 10=0. 0 
Til  =0.0 

C SELECT  WHICH  PANG  TO  USE  DEPENDING  ON  WHICH  REGION 
C YOU  ARE  IN  AND  WHICH  CASE 
IFOENPF.  EC.  I ) GOTO  230 
DA  KG  = PANG  I 

IF ( I A . GT . NR  1 ) PA  NG  = D A NG ? 

I F ( I A . GT  . ( N R UN  R ? ) ) P A KG  = D A NG  3 
IF(  I A . GT . ( N R 1 + NR  P+N  R 3 > )DA  NG  = DA  KG') 

GOTO  PRO 
230  DA  NG  = DANG  1 

IF ( I A . GT . NR  1 )DA  NG  = DA  NGE 
240  CONTINUE 

C CONVERT  A NG  INTO  EQUIVALENT  BEAM  PATTERN  ANGLE 
C FROM  PRO A PS  IDE 
XsANG-BAZ 

IF  ( BAZ . GT  . ANG)  X=X+TPI 
IF(X.GT.Pl)  X = X - T P I 
CALL  BEAM(BP, ALPHA, X, XS) 

DO  210  K=1 , NR 
NUMS=NS(K  ) 

DO  220  1 = 1,  NIIMS 

C TFST  FOR  INCLUSION  OF  SEGMENT  FOR  THIS  ANGLE 
TEMP  1=0(1 ,P, K ) + G ( I ,o,K) 

IF  ( T EM  P 1 . GE  . TP  I )T  EM  P 1 =T  EM  P 1-T  PI 
ITEST  = 1 NCLU(G ( I , F , K ) , TEN  P 1 , A NG ) 

IF (I  TEST. EQ. 0 ) GOTO  220 
C COMPUTE  MULTIPLIERS 

TEMP2=GU  , 1 1 ,K>«S(  1 , 1 ) + G(T  , IP,  K)*S(?,  1 )4C.(I  , 1?,K)»-S(3,  D 
TEMPS  = 0(7 , 1 1,K)*(S(1  , 1 )«S(1  , 1 )+S(1 , P)  ) 

1 *G  ( I , 1 ? , K ) ■ ( S ( ? , 1 ) *S  ( P , 1 ) + S ( P , 2 ) ) 

2 +G(I  , 1?,K)«(S(3,  1)«S(3,  D+S(3f  ?)) 
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C COMPUTE  ANCLE  DEVIATION  D FROM  LEFT  END  OF  DECK  FNT 
D=ANG-G(I ,e, K) 

IF(D.LT.O.O)DsD+TPI 
T1 3=5  IN (D ) 

T 1 ?=COS (D ) 

C COMPUTE  INITIAL  RANGE  ANGLE  RO  FOR  THE  GIVEN 
C ANGLE  D AND  INITIAL  GIO 
IRFLAGr  1 
0 = 0.0 

TEM. P5  = (G  ( I , 7 , K )-D)  *0.  5 
TKMPC= (0(1 ,7, K)+D) *0.5 
T EMP7  =G  ( I , 6 , K ) *S IN (T EM P5 ) /S I N (T  EM  P6  ) 

TEM  P8=G ( I , 6 , K ) * COS (T  EM  P5 ) /COS ( T EM  P6  ) 

R 0= A TA  N ( T EM  P7  ) + A TA  N ( T EM  P ? ) 

G 10  = A TA  N ( T EM  PP  ) -A  TA  N ( T F.M  P7  ) 

WI  = 0. 5*WIDTH(GI0,G(I ,2, K) ,G(I , 3, K) ,G(I , 1 , K)  ) 
IF(WI.LT. DELR)  I RFLA G=0 
R = R 0 

TEMP5=0.  0 

C LOOP  ON  RANGE  R FOR  THIS  ANGLE  AND  SEGMENT 
TEMP6=0. 0 
G J =GIO 

300  TEMP7  = 7.(R,  ANG)*BP 

TEMP8  =X J (R , Q)#PETA(Q,WJ) 

TEMP5=TEKP5+TEMPP*TEKP7 
T EM  P 6 =T  EM  P 6 +T  EM  P P * T EM  P 7 * T EM  P7 
IF ( I RFLA G . EO. 0 ) GOTO  ^50 
C INCREMENT  RANGE 

IFURFLAG.EQ.  1 )R  = R + DELR 
IF(IRFLAG.EQ.-1 )R=R-DELR 
C COMPUTE  NEW  GT,  QI , WT 

TEMP7  =G (I ,5,K)*C0S(R)  + G(I ,H , K) «S JN(R)*T  12 
TEMP8=A  COS (T  EM P7 ) 

TEM P7=A SIN (S IN (R)*T 12/SIN (TEM P8 ) ) 

G I = A TA  N ( COS ( T EM  P7 - G ( I ,7, K) ) *T A N (TEM PP ) ) 

Q=A  SIN (S IN (TEM  P7-G( 1 , 7 , K ) ) *S IN (TEM  PE ) ) 

WI  = 0. 5 fWI DTH  CGI , G( 1 , 2,  K ) , G(I  , 3 , K ) , C(  1 , 1 , K)  ) 

C TEST  TO  SEE  IF  NEW  VALUE  OF  01  PUTS  YOU  OUT  OF 
C THE  ROUTE  SEGMENT  ENVELOPE. 

TEMP8  =A  BS (0 ) -W I 

IFCTEMP8.  GE. 0. 0. AND. IRFLAG. EO. 1 ) GOTO  370 
IFCTEMP8.GE.0.0. AND. IRFLAG. E0.-1 ) GOTO  350 
GOTO  300 
370  R=RO 
Q=0 . 0 
G I =G  10 
IRFLAGr-  1 

W 1=0 . 5*WI DTH (G I , G( I , 2, K),G(I,3,K),C(T,1,K)) 

GOTO  300 

350  T 10=T 10+TEMP5*TEM P?*DELR 
T 1 1 =T  1 1+TEMP6*TEMP3*DELR 
220  CONTINUE 
C END  OF  SEGMENT  LOOP 
210  CONTINUE 
C END  OF  ROUTE  LOOP 
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XMEAK  = XMEAN-*T  10*PAIIG 
XVAH  = XVABtT  1 1'PANG 


akc=ang+dang 

I F ( A K G . GE  . T P I)  A N G = A K G - T P 1 
00  CONTINUE 

END  OF  ANGLE  LOOP 
XK  FA  Nr  I?  110  * X M FAN 
XV  A K = lt  HO  * X V A I? 

C#k  * »#  **  * ***(.;  fip  OF  MEAN  AP’D  VAPTAMCF  COM  PUT  A TTOKSr  * * * * ' * 
C»**>Ht***,**’*****itM(irk«><*tnc»*i.K  »«*#*>■*♦<»'**»***'»***#«  *♦.*»•* 

C MAIN  COMPUTATION 

C COMPUTE  REAL  PART  OF  CHARACTERISTIC  FUNCTION 

C*******!;** 

C RESET  ANGLES  AND  INITIALIZE  F PARAMETERS 
A KG: EDGE  1 
F=0 . 0 
DELFsO.Ol 
NF  = 1 01 


C * * * * « « * >:  « . > 

C LOOT  ON  ANGLE  A KG 
DO  400  I A - 1 , NT 

C INITIALIZE  SUM  ARRAY  TO  ZERO 
DO  401  KK=1,NF 
S UM ( 1 , KK)sO. 0 
SUM( ?, KK ) = 0 . 0 
401  CONTINUE 

C SELECT  WHICH  DANG  AND  WHICH  CASE 
IFUEKDF.EO.  1)  GOTO  410 
DA  KG  = DA  KG  1 

I F ( I A . GT.  NR  1 ) DA  NC>  = DA  NG  2 
I F (I  A . C,T  . ( N R UN  R 2 ) ) D A NG  = DA  NG  ? 

I F ( I A . GT . ( N R 1 +N R 2+ NR ? ) ) DA  NG - DA NG 4 
GOTO  420 
410  DA  NG  = DA  KG  1 

IF(IA.GT.NRI)  DA  KG  = DA  KG F 
420  CONTINUE 

C CONVERT  A KG  INTO  EQUIVALENT  P FA  M PATTERN  ANGLE.  X 
C FROM  PROA  PS  IDF  AND  COMPUTE  HP 
XrANG-PAZ 

IF(PAZ. GT. A KG >X  = X+TP1 
I F ( X . GT . PI )X=X-TPI 
CAI  1.  BEAM  (IP,  ALPHA,  X,  NS) 

C LOOP  ON  ROUTES 

DO  4 TO  K = 1 , NR 
NUKS=NS(K ) 

C LOOP  ON  SEGMENTS 

PO  440  1=1, RUMS 

C TEST  FOR  INCLUSION  OF  A NG  IN  TUTS  SEGMENT 
TF.KP1=GU  ,8,  K)  + G(T  ,o,K) 

IF (TEN  PI. GF. TPI )TEM P 1 =TFM ri-TPI 
I T E ST? I NC  LU (G ( 1 , 8 , K ) , TEM P 1 , ANG ) 

IF ( I TEST. EC. 0 )G OTO  440 

C COM  PUT!'  ANGLE  PE  VI  AT  ION  P FROM  LEFT  END  OF  SEGMENT 
P= A NG -G ( 1 ,r,  K) 

I F ( D . l.T . 0 . 0 ) D = D +T  P I 
T 1 3=S IN (D ) 

T 12=C0S (D ) 
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C COMPUTE  INITIAL  RANGE  RO  A l.'D  GTO 
I RF1.A  G=  1 
0 = 0.0 

TEMP5=(C.(I  , 7,K)-D)*0.5 
TEMP6=  (GO  ,7,K)+D)*0.5 
T EMP7=G  (T,6,K)*ST  N ( T ! ’ ' P5  ) /?,  1 1.'  ( T EM  P 6 ) 

T K V,  P P = G C I , 6 , K ) f C OS  (T IP'  P5  ) / CC  S (TH'  P 6 ) 

H 0= A TA  N ( T If.  P7  )-*  ATAIJ  ( T YV  P P ) 

G 1 0 = A TA  N ( T YV  PP  ) - A TA  It  ( T I-  M P7  ) 

WI  = 0. 5 * W 1 1 TII(G10,G(  T , ? , K 1 , G ( T , ' , K ) , G ( I , 1 , K ) ) 

JF(WJ . LT. CELR )IRFLAGrO 
R = R 0 
G1 =G  10 

c loop  or:  rage  h fcr  fixed  angle 

<150  TEFP7  = BP*Z(RI  AN’G) 

TEMPP  = XJ(R,  Q)  *1  ETA(Q,  WI) 

F=0 . 0 

C LOOP  ON  CHARACTERISTIC  FUNCTION  ARGUMENT  F 
DO  '160  I W = 1 , NF 
V=F*TEM  P7 

CALL  GAMA (V , CL , B1L , CEL , CIL) 

CALL  GA M A ( V , CS  , BIS,  Cl’S  , CIS  ) 

CALL  G A M A ( V ,CF,P1F,CFF,CIF) 

S U H ( 1 , I W ) = S U M ( 1 , IW)+TEKPP* (G(I , 1 1 , K)*CPL 
1 +G ( I , 1 ? , K ) *C  RS+G ( I , 1? , K ) *CPF ) 
suk(  ? , iw)  = su::(2,  iw)+tempp*(g  a , i i,k)*ctl 
1 -4 G ( I , 1 ? , K)  - CIS+Gd  , 1 ?.,K)  *CIF) 

FsF+DELF 
4 60  CONTINUE 

C INCRFMFKT  RANGE 

IF(T  RFLAG.  EC.  C)  GOTO  '170 
IF  ( I RFLAG  . 1C.  1 ) R = R + DF.LR 
I F ( I R F L A G . E 0 . - 1 ) R = R - 1 > E L R 
C COMPUTE  NEXT  Cl,  QI  AND  VJ 1 

TEMP 7 =G ( I , 5 , K ) * CCS ( R ) +G (I ,4,K)*SIN(R)*T1? 

TFKPP=ACOS  (T  EM  F‘7  ) 

T EM  T7  = A SIN  (S  IN  ( R ) *T  1 VS  I M (T  EM  PP  ) ) 

G I = A TA  N ( COS ( T EM  P7 — G ( i ,7,  K)  ) {:T  A N (TEM  PP  ) ) 

Q = A S I N ( S I N ( T EM  P7  - G ( I , 7 , K ) ) * S I N ( T EM  PP  ) ) 
v:i=0.  SMvIDTIKGT  ,G(I  , ?,  K)  ,C,(I  , ?,  K)  , G(I  , I ,K)  ) 

C TEST  TO  SEE  IF  NEW  VALUE  OF  01  PUTS  YOU  CUT  CF  THE  ROUTE 
C SEGMENT  ENVELOPE 
TFMPP  = A BS (Q ) -W  I 

IFCTEMPP. GE. 0. 0. AND. IRFLAG. EC. 1)  GOTO  4P0 
I F ( T EM  PP . GE.  0.  0.  ALT.  IRFLAG.  FC.  -1  ) GOTO  4 70 
GOTO  450 
480  R = R 0 

0=0.0 

GIrGIO 

IRFLACr-1 

W 1=0 . 5 ‘ V.’  I PTH  (G I , G(  I , 2,  K ) , G(  I , 3 , K)  , G(  1 , 1 , K ) ) 

GOTO  450 

470  T 1 4 =DELR* DA NG 

DO  4 °0  T W = 1 , NF 

PHI  ( 1 , IW)  = PI!J(1,  IW)  +S UM(  1 , TW)*T14 
HI]  (?,  IVOsPlII  (.7,  TVO+SUMC?,  IW)  *T  1 4 
490  CONTINUE 

440  CONTINUE 

C* «***E ND  OF  SEGMENT  LOOP***** 

n *5  0 
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<130  CONTINUE 

c* ****end  ok  hoiitk  loop*'**** 

ANG  =A  NG  + DA  IIC 

IF(  ANG.GE.TP1  ) ANGrA  NG-TPI 
'100  CO  NT  1 N UK 

C*«***ENP  OK  ANGLE  LOOP*  * * * * 

C COM  PUT  K RIAL  PART  OK  CIIA  If  ALTER  1ST]  C KUNCTION  KROM 
C SECOND  CHARACTERISTIC  KUNCTION 
K=0 . 0 

DO  500  JK  = 1,MK 

PUT  ( 1 , JK)  = PKI.7.»  K XT  (R IIO*  Pill  ( 1 , JK)  ) * COS  (RIIOM’H  1 O’,  .IK)  ) 

pin  (;*,  jk)=f 

K=  IT  DK1.K 
500  CONTINUE 

C 

C**''**ENP  OK  CHARACTERISTIC  KUNCTION  COMPUTATIONS***** 

C 

C OUTPUT  R K SUETS  TO  KILE  PH  I X 
C 

C WRITE  HEADER  ON  KILE  I'll  IX 
WRITE (6,  MO) 

510  KORMATdX.PfiHRESULTS  FROM  PROGRAM  CKUNK) 

WRITE  ( (',  B?0  ) 1 DN UK  , KRTC,  DELFQ 
5?0  FORMAT!  IX,  1 v’KIO.'l) 

C 


c 

c 

CONVERTS  DSTER  AND  BW IDE  TO  DEGREES  K 

ROM  RADIANS 

BSTKRr 1 PO . * DSTKR/P 
BW! DE- 1 PO. * BW 1 IT/P 

I 

I 

WRIT E ( 6 , 5 30) ESTER, 

BW  I DE 

030 

r 

FORMAT!  1 X,  ?K  IS.  '1  ) 

V 

W R 1 T E ( f> , 5'I0)TSIU  P, 

DEI. 7.  , XM KAN,  XV A R 

5 '10 

FORMAT!  1 X , K 1 1 . '1 , T 
WRITE  (0,  S'lMNK  , PM 

PO.  1 '1 ) 

K 

5 '15 

FORMAT! IX,  IN, K ?0. 
WRIT!  (0,  5 50  ) (I'll  1 ( 1 

I'D 

, K)  , Pll  1 (?,  K)  ,Kr.1 

, NK  ) 

550 

FORMAT! IX, EDO.  IN,  3 

X , EDO.  I '1  ) 

c 

c 

C REWIND  I'll. IS 
REWIND  I 
REWIND  P 
REWIND  3 
REWIND  (> 

C 

STOP 

END 
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3.5  Subprograms 

This  section  describes  the  subroutine  and  function 
subprograms  that  are  called  from  the  main  program  of  CFUNK. 

3.5.1  Function  Z 

This  function  subprogram  is  used  to  return  a value  of 
transmission  loss  to  the  main  program  from  input  values,  of 
range  and  azimuth.  The  array  TL  containing  the  transmission 
loss  table  is  stored  and  labeled  common  XLOSS.  The  calling 
sequence  is: 

Function  Z (R,  D) 

■where  Z = returned  value  of  transmission  loss  in  power 

R = range  as  earth  centered  angle  in  radians. 

which  is  converted  to  n.in.  in  this  function 

D = azimuth  in  radians. 

Figure  8 presents  a flowchart  for  Function  Z.  Upon 
entry,  five  delimiting  angles  A1  to  A9  in  radians  are  sot 
to  divide  up  the  azimuth  into  the  five  regions  where  each 
of  the  columns  of  array  TL  are  to  be  used.  For  example, 
if  the  transmission  loss  Is  given  at  the  five  azimuthal 
angles  69  , 132°,  285°,  294°  and  391  (where  69  is  the  first 
column  of  TL  and  39.1°  the  last)  then  choose  A1  = 30*'  , A2  = 100. 9'  , 
A3  = 208.9  , A4  = 289. 91  and  A9  = 322.5°  as  the  delimiting 
angles.  Note  that  the  programmer  must  set  the  values  of  Al 
to  A9  in  the  program  before  execution. 

Next,  the  function  determines  the  row  index  1 from 

I = (range  in  n.in.)/10  + 1 

The  function  sets  Z = 0 and  tests  to  see  if  the  index  1 > 390; 
that  is  if  the  range  is  greater  than  3900  miles.  If  this  is 
the  case,  the  function  returns  to  the  main  program. 

If  not,  tiie  program  next  tests  to  determine  the  two 
values  of  delimiting  angles  that  the  angle  P ties  between, 
and  thus  selects  the  appropriate  column  index  J . Once  l 
and  J are  found,  the  function  then  computes  a value  of 
transmission  loss  Z by  linear  interpolation  on  the  range. 

The  value  of  Z is  then  returned  to  the  main  program. 
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FIGURE  8 Flowchart  Function  Z (R,D). 
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The  listing  of  the  FORTH  AN  code  l'ot*  Fund  Ion  Z is  as 
follows : 


FUNCTION  7.  ( R , D) 

COMMON/XLOSS/TL  ( 3 r>  1 , r> ) 

C 

C RETURNS  A TRANSMISSION  LOSS  FROM  THE  TABLE  TL 

C TL  IS  THE  INPUT  TRANSMISSION  LOSS  TABLE  FOR  FIVE 
C AZIMUTH  ANGLES.  TEN  N.M.  INCREMENTS  IN  RANGE 
C 7 = RETURNEE  VALUE.  OF  TRANSMISSION  LOSS 
C R = RANGE 

C D = ANGI  E FROM  NORTH 
C 

RAPS(V)=V*  ( .017H5329?ri) 

R HO =34 37. 74677 1 

C CONVERT  RANGE  FROM  RADIANS  TO  N.M. 

RX  = R*R  HO 

C SET  DELIMITING  ANGLES 
A I = R A DS  ( 3 0 . ) 

A?  = R APS ( 1 00.5  ) 

A 3=R  APS (?  OR . 5 ) 

A •)  = R A PS  ( ? 89 . 5 ) 

A5=RAPS(3?2.5) 

C 

C FIND  RANGE  INDEX  I (ROW) 

I =1  NT ( R X/ 1 0. ) + l 
Z=0. 

IFU.GE.350)  RETURN 
C 

C FIND  ANGLE  INDEX  J (COLUMN) 

IF  ( D . LT.  A.?.  AND.  P.  GE  . A 1 ) J = 1 
IF(D.LT. A3. ANP.D.CE. A?)  Jr? 

IF ( D . LT. AH. ANP.P.GE. A ' ) Jr3 
IF( D . LT. AS. ANP.P.GE. AH)  J=H 
IFCD.LT. A1.0R.P.GT. A5)  J =‘5 
C 

C INTERPOLATE  OVER  RANGE 
C 

SLOPE=(TL(I +1 ,J)-TL(I ,J))/10. 

RI= ( FLOAT ( I )-I . ) *10. 

Z=TL(I  , J)  + (RX-R  I ) * SLOPE 

RETURN 

END 
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3.5.2  Function  BETA 

The  computation  of  the  characteristic  function  4>y(w) 
requires  a probability  density  function  for  the  ship  position 
across  the  route.  In  Section  2 of  BBN  Report  No.  3 1> this 
density  function  is  denoted  for  tin  ith  route  by  1’  (q),  where 

q is  the  earth  cent.ei'ed  angle  in  radians  ropresent-i 
lng  across  route  s.hip  position  measured  from  the  center  of 
the  route. 

Define  the  width  of  the  route,  at  a given  point  along 
the  route,  as  W where  W is  an  earth  centered  angle  in 
radians.  Define  the  half-width  of  the  route  as 


qo  = W/R 

Therefore,  the  probTrbility  of  a ship  being,  positioned  outside 
the  interval  -q  < q < + q is  zero. 

Since  exact  information  of  ship  distributions  across  the 
world's  shipping  routes  is  not  available  at  present,  a probability 
density  function  must  be  chosen  to  model  the  uncertainty  in 
ship  position  within  the  route  width  W.  An  appropriate  choice 
is  felt  to  be  the  BETA  probability  density  function  with  the 
parameter  values  both  set  equal  to  two.  This  density  function 
is  symmetric  about  its  mean  value,  which  if  chosen  as  the 
center  of  the  route,  does  not  favor  one  side  of  a route  over 
the  other . The  density  function  to  be  used  for  ships'  posit  Ion 
across  the  route  is 


Figure  9 shows  a plot,  of  the  density  function. 
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Function  BETA  was  written  in  FORTRAN  to  return  a value 
of  (q)  from  inputs  q and  q . The  calling  sequence  for 
yi  this  function  is  BETA  (Q,  QO)  where 


BETA  - f (q);  returned  value  of  density  function. 
yi 

Q - q in  radians 


QO  - 1/2  route  width  in  radians. 


Figure  10  presents  the  flowchart  for  Function  BETA.  The 
function  listing  is  as  follows: 


FUNCTION  BETA(O.OO) 

C TRANSVERSE  SHIP  POSITION  DENSITY  FUNCTION 
C BETA  = RETURNED  VALUE  OF  DENSITY  FUNCTION 
C Q r ACROSS-ROUTE  EARTH  CENTERED  ANGLE  (RADS.) 
C QO  = ONE  HALF  THE  ROUTE  WIDTH 

C 

BETA=0 . 0 

IF(Q.LT.-QO.OR.Q.CT.QO)  RFTURN 
BETA= 1 . 0 

IF(QO.EQ.O.O)  RETURN 

X =Q/Q0 

X=1 .0-(X«X) 

BETA=(1.875*X*X)/(2.*Q0) 

RETURN 

END 
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FIGURE  10  Flowchart  for  Function  BETA 
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3.5.3  Function  INCLU 

This  function  is  used  to  determine  if  an  ancle  C is 
included  between  ancles  A and  B.  Ancles  A,  B and  C are 
azimuthal  angles  in  radians.  It  is  assumed  that  A is  always 
counter  clockwise  from  B.  The  calling  sequence  is 

INCLU  ( A ,B , C ) 

where 

the  returned  value  of  INCLU  = +1  if  C is  included 
between  A and  B,  and  zero  otherwise. 

Figure  11  presents  a flowchart  for  FUNCTION  INCLU.  The  listing 
of  the  function  is  as  follows: 


FUNCTION  I NCLU ( A , B, C) 

C FUNCTION  DETERMINES  IF  ANGLE  C IS  INCLUDED 
C BETWEEN  ANGLES  A AND  B. 

C INCLUr  0 IF  NOT  INCLUDED 
C INCLU=  1 IF  INCLUDED 
C 

INCLU=0 

IF(B.GT. A. AND. C.CE.A.AND. B.GE.C)  INCLU=1 

IFCA.GT.  B.  AND.  B.GE.C)  INCLIJ=1 

IF(A.GT. B. AND. C.GE.A)  INCLU=1 

RETURN 

END 
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3.5.4  Subroutine  GAMA 

This  subroutine  computes  the  real  and  imaginary  paid  or  the 
characteristic  function  of  the  gamma  probability  density 
function.  The  gamma  probability  density  function  is  used  to 
statistically  model  the  radiated  noise  from  merchant  ship-'. 

The  characteristic  function  is  given  by 


4>  ( to  ) 


1 


where 

C = M/a2 
b + 1 = M2/o2 
M = Mean 


o2  = Variance 


and  j = iTIT 


The  two  parameters  C and  b + 1,  computed  in  the  main  program 
from  the  source  mean  and  variance,  are  inputs  to  the  sub  rout ine  . 
The  characteristic  function  can  be  written  in  terms-  of  its- 


real 

and  imaginar 

y part  s as 

Real  <Kw)  = 

C 

b + 1 | 

COS 

(b  t 

1) 

Arc  tan  {o>  C}\ 

c2  + CO2 

( 

and 

r 

IM  4>U)  = 

c 

b + 1 

sin 

1 (b  + 

1) 

Arc  tan  (o>  C ){ 

_V  C2  + uf 

| 
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The  cal. 1.1  nr,  sequence  for  the  subroutine  is  GAMA  (W,  C , 

D,  CK,  Cl)  where  the  inputs  are: 

W = argument  u v f characteristic  funct  ion 

C = parameter  C of  the  gamma  probability  density  t’unct  Ion 
D = b + 1 parameter 
and  the  outputs  are: 

CK  = real  part  of  characteristic  function 

Cl  = imaginary  part  of  characteristic  function. 

p i purr  i presents  flowchart  for  Subroutine 
listing  of  the  program  is  as  follows: 


SUPROUTINF  OA M A ( W ,C, D, CR, Cl) 

C COMPUTES  THE  CHARACTERISTIC  FUNCTION  OF  GAMMA  TENSITY 
C Wr  FREQ.  IN  RAPS. 

C C=  PARAMETER  C OF  DENSITY 
C D=  P+ 1 PARAMETER  OF  DENSITY 
C CR=  R FA I.  PART 
C CI=  IMAGINARY  PART 
XrW/C 

A = ( 1 . O/SOR T ( 1 .0  + X*  . ) ) * * D 
Y = D * A TA  N ( X ) 

C R=A*  COS ( Y ) 

C I =A*  S IN  ( Y ) 

RETURN 

END 
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3.5.5  Subroutine  Beam 


In  Section  4.1  of  BBN  Report  3653  a normalized  beam 
pattern  model  was  described  in  detail.  The  normalized 
beam  pattern  values  is  needed  in  the  main  program,  w here 
it  is  used  to  weight  the  transmission  loss  values. 


This  subroutine  was  written  to  compute  the  normalized 
beam  pattern  model.  The  calling  sequence  is  SUBROUTINE 
BEAM  (BP,  ALPHA,  X,  XS)  where 

BP  = value  of  normalized  beam  pattern  to  be  returned. 


ALPHA  = value  of  frequency  dependent  constant  which  is  a 
program  input . 

X = angle  in  radians  with  respect  to  array  broadside 
in  radians,  which  is  an  input  to  the  main  program. 


Figure  13  presents  the  flowchart  for  Subroutine  BEAM.  A 
listing  is  as  follows: 


C 

SUBROUTINE  BEAM( BP, ALPHA , X, XS) 

PI =3 . 14159265 

B ETA  = A L PH  A * ( S IN  (X  ) -S  IN  ( X S ) ) 

IF C BETA. EO. 0. ) GO  TO  150 
A VAL=ABS(  BETA ) 

IF  ( AVAL  . GT . 2.  ) GO  TO  100 
IF ( AVAL . EQ . 1.)  GO  TO  200 

TEMP  1 =S  IN  (P I *PETA  ) 

T EM  P 1 =T  EM  P 1 * T EM  P 1 

TEMP2= 1 . 0/( (PI*BETA ) *( 1 . 0-B ETA* B ETA ) ) 
BP=T  EM  P 1 *T  EM  P2*T  EM  P2 
GO  TO  250 

100  BP=1.0/((PI*BETA)*(1.0-BETA*BETA)) 

BP=0 . 5*  BP*  BP 
GO  TO  250 
150  BP= 1 . 

GO  TO  250 
200  B P=0 . 25 

250  RETURN 

END 
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4.  PROGRAM  DLNS 

The  final  computer  program  In  the  implementation  of  the 
Ambient  Noise  Model  is  Program  DENS . 

'I'he  Ambient  Noise  Model  produces  as  output  statistical 
measures,  of  the  noise  power  at  the  output  of  a beumformer. 

The  measures  that  are  of  interest  are  the  probability  density 
function  of  output  noise  power  (in  power),  the  density 
function  of  output  noise  power  (in  db  re/watt),  and  the 
distribution  function  of  output  noise  power  (in  db). 

Program  DENS  calculates  the  probability  density  function 
for  the  beamformer  output  noise  power  from  the  real  part  of  its 
characteristic  function.  This  inverse  transformation  w . 
derived  and  discussed  in  Section  3.1  of  I3BN  Report  3Gl  3.  One 
density  function  is  found  as  a function  of  power  (waits)  and 
another  as  a function  of  db  (re/watt).  Next  , tin  pr.  ■ • 
calculates  the  distribution  function  (in  db)  frem  the  • ■ v 

function  in  power.  The  program  reads  input  data  from  tin 
previously  created  disk  file  P111X  and  prints  the  result, 
at  the  terminal. 

4.1  Input  Data  File 

The  input  data  for  Program  DENS  is  read  from  disk  file 
PHIX.  The  contents  and  structure  of  the  logical  records,  in 
file  PllIX  were  described  in  Section  3.3  above. 

4.2  Main  Program 

The  main  program  of  DENS  performs  the  following  functions: 

1)  initialize  arrays  and  set  program  constants 

2)  reads  input  data  from  file  PHIX 

3)  computes  the  required  probability  density 
and  distribution  functions 

4)  prints  the  results  at  the  terminal. 

Program  DENS  calls  three  subroutines  from  the  main 
program  during  its  execution.  These  subroutines  are 

• IPT  - Calculates  the  inverse  fourier  transform,  CJX, 
for  a given  value  of  X from  the  real  part  of  the 
characteristic  function  in  array  R.  Array  R,  the 
number  of  sample  points  JIM,  and  the  sample  point 
spacing  DELW.1.  arc  stored  in  the  labeled  common 
area  COM . 
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• P0W2DB  - Converts  a probability  density  function 
in  power  into  a probability  density  function  in 
db/re/watt . 

• DIET  - Calculates  the  distribution  function  in 
power  from  an  input  density  function  in  power. 

These  subroutines  are  described  in  Section  4.5  below. 

4.2.1  Common  Areas 

Program  DENS  makes  use  of  one  common  area  called  COM. 

COM  contains  the  real  array  R,  a real  constant  DELW1,  and 
an  integer  constant  JIM.  COM  is  used  by  the  utility  sub- 
routine IFT.  Array  R contains  the  real  part  of  a character- 
istic function.  JIM  is  the  number  of  sample  points  in  the 
characteristic  function,  and  DELW1  is  the  sample  point  spacing. 

4.2.2  Description  and  Flowchart 

Figure  14  presents  a flowchart  for  Program  DENS.  The 
program  starts  by  initializing  constants  and  arrays.  Next, 

DENS  reads  data  from  the  incut  data  file  PHIX.  These  items 
are  the  header;  an  identification  number,  center  frequency, 
and  bandwidth;  the  beam  steering  angle  and  major  lobe  width; 
the  total  mean  number  of  ships,  the  probability  of  no  noise, 
the  mean  noise  power,  and  the  variance  in  noise  power;  the 
number  of  sample  points  and  the  sample  spacing  of  the  character- 
istic function;  and  finally,  the  real  part  of  the  characteristic 
function  stored  in  array  R. 

Now,  the  first  of  three  functions  is  calculated.  The 
first  function  is  the  probability  density,  f(X),  for  the  noise 
power  at  the  output  of  a beamformer  (in  power) . This  func- 
tion starts  at  X = 0 . The  calculation  proceeds  in  two 
steps;  first,  the  contribution  due  to  an  impulse  at  zero 
frequency  is  removed  by  subtracting  DELZ  from  each  element  of 
the  characteristic  function,  R(J).  Second,  after  entering  a 
DO  loop  on  X,  the  result  from  Step  1 is  inverse  fourier 
transformed  using  utility  subroutine  IFT.  IFT  produces  a 
value  for  f(X)  for  each  passage  through  the  loop.  In  addition, 

X is  incremented  by  DELX  on  each  passage.  The  values  of  f(X) 
and  X are  stored  in  FX(l,J)  and  FX(2,J),  respectively.  The 
Index  J runs  from  1 to  NX. 


I 
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FIGURE  1 


4 Program  DENS  Flowchart. 
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The  second  function  that  is  calculated  by  DENS  is  the 
probability  density  in  db  (re/watt)  of  the  noise  power  at 
the  beamformer  output.  It  is  found  by  using  the  utility 
subroutine  P0W2DB  to  convert  the  density  function  in  power, 
f(X),  into  a density  function  in  db , f(y).  However,  the  point 
X = 0 watt  is  skipped;  therefore,  f(y)  has  one  less  sample 
point  than  f (X)  . The  values  of  f(y)  at  the  sample  points 
y are  stored  in  FY(l,J),  while  the  sample  points  are  stored 
in  FY(2,J) . 

The  third  function  that  is  calculated  by  DENS  is  the 
distribution  function  (in  db  re/watt)  for  noise  power  at  the 
beamformer  output.  This  function  is  found  by  using  the  utility 
subroutine  D1ST,  which  produces  a distribution  function  in 
power  from  a density  function  in  power.  The  distribution  func- 
tion (in  db)  is  calculated  in  two  steps;  first,  DIST  is  called 
to  operate  on  FX . This  produces  DF(I,J),  where  the  DF(1,J) 
are  the  sampled  distribution  function  (in  power),  and  the 
DF(2,J)  are  the  sample  points  Sj . Second,  the  abscissa  of 
the  distribution,  DF(2,J)  are  converted  to  db  re/watt,  and 
stored  back  in  DF(2,J)  (note  the  point  X = 0 is  skipped). 

Output  data  are  printed  at  the  terminal.  This  includes 
descriptive  information  (such  as  center  frequency,  band  width, 
beam  steer  angle,  etc.)  that  is  pertinent  to  the  case  at  hand, 
followed  by  three  tables  of  the  density  and  distribution 
functions . 

The  maximum  number  of  sample  points  is  fixed  in  the  program 
as  301  for  R and  for  FX  (i.e.,  NX  = 301).  While  the  maximum 
number  for  FY,  and  DF  are  one  less,  300.  The  sample  point 
spacing  of  f(X),  DELX  is  equal  to  101*  . These  limitations  can 
be  changed  easily  by  altering  COMMON,  DIMENSION,  and  DATA 
statements . 

4.3  Output 

The  output  from  Program  DENS  is  printed  at  the  computer 
terminal.  This  is  the  final  output  of  the  Ambient  Noise  Model 
for  the  particular  case  of  interest . 

The  program  prints  the  identification  number  for  the 
particular  route  and  sensor  geometry;  the  frequency  band  and 
the  center  frequency;  the  beam  steering  angle  and  beam 
width  of  the  main  lobe  of  the  array;  the  total  number  of 
ships;  the  probability  of  no  noise;  the  mean,  variance  and 
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standard  deviation  of  benmformor  output  noire  power  in  units 
of  watts;  the  probability  density  function  of  noise  power 
in  watts;  the  probability  density  function  of  noise  power 
in  db ; and  the  probability  distribut  ion  function  of  noise 
power  in  db.  Table  o shows  the  structure  of  the  printed 
output  from  Program  DENE. 

4.4  Glossary  and  Program  Listing 


The  fo .1  lowing  glossary  contains  definitions  of  the  FORTRAN 
constants  and  variable  names  used  in  Program  DENS.  The  names 
are  presented  in  alphabetical  order. 


BSTER 
BW  l DK 


DELE 


DF.LFQ 


DKbW.1 


beam  steer  angle  (dog.) 
width  of  major  lobe  (deg.) 

sample  point  spacing,  'a  oharacterlst lc  function 
(rad i ans/sec ) . 

beamformer  bandwidth  (Its). 

sample  point  spacing;  used  in  subroutine  I FT 
( rad i ans/sec ) . 


DELX 

DELZ 

m-'(i,J) 


DF ( 1 , J ) 


FREQ 


FX  (.1,3) 

FX(2,.l) 
FY (1,3) 


FY(2 , J ) 


sample  point  spacing  of  density  function  1'tX). 
probability  of  no  noise. 

distribution  function  of  beamformer  output  noise 
at  sample  point  J. 

sample  point  J of  distribution  (db/re  watt). 

beamformer  center  frequency  (lls). 

probability  density  (in  power)  of  benmformor 
output  noise  at  point  X; 

sample  point  X;  in  probability  density  (watts) 

probability  density  (in  db  re/ watt)  of  beamformer 
output,  noise  at.  point  y . 

sample  point-  y;  in  probability  density  (db  re/watt) 
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TABLE  6 STRUCTURE  OF  THE  OUTPUT  OF 
PROGRAM  DENS 

(NOTE:  Each  line  represents  one  line  of  printout) 


10  NUMBER  - 

FRh'QUf  NO  - HZ 

CENTER  FREQUENCY  = HZ 

Bt  AM  STEERING  ANGLE  - DI  G 

BEAM  WIDTH  (MAJOR  LOB  I Y = DEG 

TOTAI  AVG.  NUMBER  Ol  SHIPS  = _ 

PROBAB 1 1 I TV  OF  NO  NOISE  =_ 

MIAN  BLAMFORMER  OUTPUT  NOISI  POWER  = WATTS 

VARIANCE  OF  BLAMf'ORMER  OUTPUT  NOISE  POWER 
WATTS*  * P 

STAND.  01 V.  IN  OUTPUT  NOISE  POWER  = 

WATTS 


NX-1 
l INFS 


NX-1 

LINES 


NX  ( 
LINES 


PROBAB 11  1 TY 

POW1  R (WATTS) 

OLNS 1 TY 
(IX) 

(X) 

i 

PROBAB 1 1 1 TV 

POWER  (OB) 

OENS 1 TY 

(Y) 

(FY) 

(RE/WATT) 

OISTRIBUT ION 

POWER  (OB) 

FUNCT ION 

(Y) 

(OF) 

(RE/WATT) 

‘ 1-0 
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IDNUM 

- 

Identification  number  for  thin  cane. 

J 1M 

- 

number  of  sample  points;  used  In  subroutine 

I FT  = NF. 

JN 

- 

number  of  sample  points  -1;  used  Internally  to 
skip  X = 0 in  tlb  conversions. 

NF 

- 

number  of  sample  points  in  character! st tc  function 

NX 

- 

number  of  sample  points  In  f(X). 

PHIX 

- 

input  dat a file. 

R(J ) 

- 

real  part  of  characteristic  function;  used  in 
subroutine  1 FT . 

TUMP ( 1.3) 

- 

working  array. 

TSMI P 

- 

total  mean  number  of  ships. 

XMKAN 

- 

mean  beamformcr  output  noise  power  (watts). 

XSTD 

- 

standard  deviation  of  beamformcr  output  noise 
power  (watts). 

XV  AH 

- 

variance  of  beamformcr  output  noise  power  (watts2) 

The  listing  of’  Program  DENS  is  at?  follow;?: 


PROGRAM  PENS  (OUTPUT,  PHIX,  TA  PF  1=PHIX) 



C THIS  PROGRAM  CALCULATES  THREE  FUNCTIONS  THAT  ARE  IJSFP  IN  THE 
C AMBIENT  NOISE  MOPEl.  . FIRST,  IT  FINDS  THE  PROBABILITY  DENSITY 
C FUNCTION  ( IN  POWER  > FROM  ITS  CHARACTERISTIC  FUNCTION.  SECOND, 
C THIS  DENSITY  FUNCTION  IN  POWER  IS  CONVERTED  INTO  A DENSITY 
C FUNCTION  IN  DD/  RE  1 WATT.  THIRD,  THE  DISTRIBUTION  FUNCTION 
C ( IN  DB  ) IS  CALCULATED.  WRITTEN  21 JULY  1978  BY  W. SCOTT. 

C 

COMMON /COM/ R ( 3 0 1 ) , DELW1 , J 1M 

DIMENSION  FX(2, 301) , FY(2, 700 ) , DEC  2, 701) , TEMPI 2, 7 00) 

DATA  R/7  01  *0. 0/  , NX/7  01  / , PELX/1 . 0E0‘1/  , TF.M P/6 00*0. 0/ 

DATA  FX/602*0. 0/ , EY/600*0. 0/ 

REWIND  1 

C READ  IN  TRASH 

REAP(0 1 , 100) 

100  FORMAT( 1 X, 26H 


) 
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C READ  TDNUM,  FREQ,  DFLFO 

R FA  P( 0 1 , 150)  IDNUM,  FRFC,  DELFQ 
150  FORMAT ( 1 X, 13, ?F  10. D ) 

C READ  BEAM  ST FF HI  NO  ANGLE  AND  MAJOR  LOBE  WIDTH 
READ(01,P00)  PSTFR,  PW I DE 
POO  FORMATt 1 X,  PF  1 5.  ** ) 

C READ  TSHIP,  DEL.'!,  XMEAN,  XVAR 

READ  (01,250)  TSHIP,  PELZ , XMFAN,  XVAR 
250  FPRMAT( 1 X, F 1 1 . *( , 3EP0. 1 « ) 

C READ  NF,  AND  DFLF 

REAPCO 1,300)  NF,  PELF 
300  FORMAT!  1 X,  ID,  F, 20.  1 D ) 

C READ  IN  THE  CHARACTERISTIC  FUNCTION 
R F A D ( 0 1 , 350)  ( R (K ) , K= 1 , NF) 

350  FORMAT! IX,  EDO.  1*1,/) 

C SET  THE  NUMBER  OF  SAMPLE  POINTS  AND  THE  FREQUENCY  SPACING 
C IN  THE  CHACTERISTIC  FUNCTION. 

JIM  = NF 
DELW1  = DELF 

C REMOVE  CONTRIBUTION  FROM  DELTA  FUNCTION  AT  FREC=0.0  . 

DO  1»00  1 = 1 , NF 

R(I)  = R (I ) -DELZ 
**00  CONTINUE 

C CALCULATE  THE  PROBABILITY  DENSITY  IN  POWER 
X a 0.0 

DO  *150  J = 1 , NX 
CALL  IFT(GX,X) 

FX( 1 , J ) = GX 
FX(2,J)  = X 
X = X+DELX 
*<50  CONTINUE 


C FIND  TIIF  PROBABILITY  DENSITY  IN  DB/RE  1W  ATT.  (SKIP  X = P.P) 


JN 

DO 


= NX-1 
500  J = 
TEMP 


MP(1, J 
TEM  P(  P , J 

500  CONTINUE 


1 , JN 

J) 

.0 


= F X ( 1 , J + 1 ) 
= F X ( P , J 4 1 ) 


CALL  POW ?D B ( T EM  P , JN , F Y ) 


C FIND  THE  DISTRIBUTION  FUNCTION  IN  DB/RF  1WATT 
C FIRST  FIND  THE  DISTRIBUTION  FUNCTION  IN  POWER;  THEN, 
C CONVERT  THE  ABSCISSAE  TO  DB . 
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CALL  DISKFX,  NX,  DF  , DELZ) 

DO  550  K = 2,  NX 

DF( 2 , K ) = 1 0.*AL0G10(DF(2, K)  ) 

550  CONTINUE 

XSTD  = SQRT(XVAR) 

C PRINT  DATA  AT  TERMINAL 

PRINT  1900 
PRINT  650 

650  FORMAT ( 1 X, 4 OH  THE  FOLLOWING  DATA  ARE  THE  RESULTS  OF) 

PRINT  700 

700  FORMAT ( 1 X, 36H  PROGRAM  DENS,  WRITTEN  21JULY78  BY) 

PRINT  750 

750  F0RMAT(1 X, 1 1H  W. SCOTT.) 

PRINT  1900 
PRINT  800,  IDNUM 

800  FORMATOX,  12HID  NUMBER  = ,17) 

PRINT  850,  FREO 

850  FORMATOX,  12HFRE0UENCY  = , F10.4.3H  HZ) 

PRINT  900.DELFC 

900  FORMATC IX, 19H CENTER  FREQUENCY  = ,F10.4,3H  HZ) 

PRINT  950,  BSTER 

950  FORMATO X, 22HBEAM  STEERING  ANGLE  = , F15.4.4H  DEG) 

PRINT  975,  BWIDE 

975  FORMATOX,  19HMAJOR  LOBE  WIDTH  = ,F15.4,4H  DEG) 

PRINT  1000, TSHI P 

1000  FORMAT( IX, 29H TOTAL  AVG.  NUMBER  OF  SHIPS  = , E 1 1 . 4 ) 

PRINT  1050,  DELZ 

1050  FORMATO  X,  26HPROB ABILITY  OF  NO  NOISE  = ,E12.6) 

PRINT  1100,  XMEAN 

1100  FORMATO  X,  37HMEAN  BEAMFORMER  OUTPUT  NOISE  POWER  = , 

1 E12.6.6H  WATTS) 

PRINT  1150 

1150  FORMATOX,  HHHVARIANCE  OF  BEAMFORMER  OUTPUT  NOISE  POWER 
PRINT  1175,  XV  A R 

1175  FORMATOX, E12. 6, 9H  WATTS**2) 

PRINT  1200 

1200  FORMATO X, 35HSTAND.  DEV.  IN  OUTPUT  NOISE  POWER  =) 

PRINT  1225,  XSTD 

1225  FORMATOX,  E12. 6, 6H  WATTS,////) 

PRINT  1250 

1250  FORMATOX,  27H  PROBABILITY  . POWER  (WATTS ) ) 

PRINT  1300 

1300  FORMATOX, 22H  DENSITY  . (X)) 

PRINT  1350 

1350  FORMATOX,  14H  (FX)  .) 

PRINT  1400 

1400  FORMATO  X.27H ) 

PRINT  1500,  ( F X ( 1 ,J) , FX(2 , J ) , J = 1 , NX) 


4-9 


t 3654 


Bolt  Beranek  and  Newman  Inc. 


1500  FCRMAK1  X,  E12.6, 3»  . ,E12.6) 
PRINT  1000 
PRINT  1550 

1550  F CRM. AT (IX,  27H  PROBABILITY  . 
PRINT  1600 

1600  FORMAT! 1 X,  22H  DENSITY 
PRINT  1650 

1650  FORMAT! 1 X, 25H  (FY) 

PRINT  1 '400 

PRINT  1700,  (FY( 1 , J) , FY(2, J) 
1700  FORMAT! 1 X, E 12. 6, 3H  . rF7.2) 
PRINT  1900 
PRINT  1750 

1750  FORMAT! 1 X, 25H  DISTRIBUTION. 
PRINT  1800 

1800  FORMAT! 1 X, 22H  FUNCTION 
PRINT  1850 

1850  FORMAT! 1 X, 25H  (DF) 

PRINT  1 4100 

PRINT  1700,  ( DF! 1 , J ) , DF ( 2 , J ) 
PRINT  1900 
1900  FORMAT! 1 X,///) 

REWIND  1 

STOP 

END 


POWER (DB)  ) 
( Y ) ) 

! R E/WATT ) ) 
J=1, JN) 

POWER (DB)  ) 
(Y)  ) 

(RE/WATT)  ) 

J =2 , NX) 
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4.5  Subroutines 

This  section  describes  the  three  subroutines  that  are 
called  from  the  main  program  of  1 > K i J . . 

4.5.1  Subroutine  1FT 

In  Section  3.1  of  BBN  Report  lot  j,  a mum  r leal  rm.  the. I 
was.  developed  to  evaluate  the  special  form  of  the  inverse 
Four'cr  transform  given  by: 


f(x) 


+ 03 

O 

— R(w)  costor.  dx, 
0 


foi*  x > 0,  where  R ( *•)  is  the  real  part  of  the  eh  iracter- 
istic  function  of  the  probability  density  function  f(X)  of  the 
positive  real  random  variable  X. 


The  purpose  of  this  section  is  to  document  the  Subroutine 
1FT,  which  implements  this,  numerical  method. 


The  calling  sequence  is  Subroutine  IPT(FX.X)  where  X is  the 
value  of  the  average  square  pressure  and  PX  is,  the  value  of 
the  probability  density  function  returned  by  the  subrout  Uu  . 

Data  ls  passed  from  the  main  program  to  the  subrout ine  in 
a common  storage  area  labeled  COM,  which  contains  the 
following : 

R(301)  - an  array  containing  the  real  part  of  the 
characterist ic  fund  ion, 

DELW1  - the  integration  step  slue  in  the  w domain, 

JIM  - the  number  of  sample  points-  in  the  w domain. 


Figure  lb  presents,  a flowchart  for  SUBROUTINE  IFT. 

Upon  entry,  the  subroutine  tests  to  see  if  X is  equal  to  zero. 
If  X is  equal  to  zero  then  f ( 0 ) is  computed  by  a numerical 
approximation  to  the  equation 


f(0) 


x 


4-11 
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n 


co.vruu  »x  jkpm: 


f(o)  - :A‘'  (*n)  R<0)) 
}" 


R(n) 


(n  by  step*  of  2) 


COMPUTE  tx  from: 

1»R(2) 


u.Xj  x 
N -2 


I |[R(n-2)-4R(rt-|)^R(«)  .om-^x  - : V'  .\RJn-J  ' ’ ' ( : s ; n U1  xl 

n*‘  ' ° AtdX  ■*  *'  ) 

(n  90c\  by  steps  of  ?) 


r«(s  ?)  -iRltl-l  ).  >R(N)1 

A ' U(N  .M-.-RiS  D<R(hM1 

1 ,A»*J  i 

l \J>  J 

^ RETURN^) 


(n 
si 

?R(n*t ) 


FIGURE  15  Flowchart  of  Subroutine  1FT. 
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3- 'in 


developed  In  Section  3.1  of  BiUJ  Report  3653.  H X > 0,  then 
f (X)  Is  computed  by  the  appropriate  formula  developed  in 
Section  3.1  ol‘  BBN  Report  3653. 

A listing  of  Subroutine  11' T is  as  follow;-; 


SUBROUTINE  1FT(FX,X) 

C COMPUTES  INVERSE  FOURIER  TRANSFORM 
COMMON/COM/ R (301 ) , PELW  1 , JIM 
P I - 3 . 1 U 159265 

K 1 = JIM-2 

IF  ( X . NE.  0.  0)00  TO  ISO 
FACT1 =2. 0/(3. *PI  ) 

SUM  1=0. 

DO  100  K = 1 , K 1 , 2 

SUM 1=SUM1+R(K)+2.*R(K+1  ) 

100  CONTINUE 

SUM  1 =DEI.W  1 * ( R ( J 1 M ) - R ( 1 )+2.*SUM1) 

FX=F  ACT  1 * SUM  1 
GO  TO  300 
C 

150  C = C OS  ( 2 . * D E LW  1 * X ) 

S=SIN(2.*DELW1*X) 

Y2  = 1 . 

Y 1 =0 . 

SUM  1=0. 

F A C T2  = 1 ,/(PI*X*X) 

FACT3=F  ACT2/X 
DO  200  K = 3 , K 1 , 2 

FACTU=R(K-2)-n. *R(K-1 ) + 6.*R(K)-'l.*R(K+1 ) +R (K  +2  ' 

FACTS  = (2.*R(K-2)-U.*R(K-1  ) + R . * R (K  + 1 )-2  . * R (K+2  ) ) / ( DF.LVf  1 *X  ) 
SN=Y  1 *C+Y2*S 
CN=Y2*C-Y 1 #S 
Y1=SN 

Y 2=C  N 

SUM  1 =SUM 1+F  ACT4  *CN-F  ACT5*S  N 
200  CONTINUE 
C 

SN=Y  1 *C+Y?*S 

CN=Y2*C-Y1*S 
Y 1 =S  N 

YP=CN 

S UM 1 =S UM 1 + 3 , * R ( 1 )-'l.*R(2)+R(3) 

SUM  1 =S  UM 1 + (R ( J TM-2 )- A . * R ( J 1M- 1 ) + 3 . * R ( J 1M ) ) * CN 
SUM 1=SUM1*F ACT2/DELW 1 

TEMP=R(J 1M-2 )-2  . * R ( J 1M-1 ) + R(J 1M) 

TEM  P=2 . * ( R ( J 1M1 * X * X -T  EM P/ ( DELW  1 »DELW  1 ) ) 

FXrSUM 1 +TKM  P*S  N*  F AC  T 3 

C 

300  RETURN 
END 
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Subroutine  l’0W2l)B 

This  subroutine  converts  n probability  dens  1 1 y I'unct  Ion 
ol'  a random  variable  in  power  t o a probability  density 
I'unct  ion  of  a random  variable  in  db. 

bet  Y be  a random  variable  in  units  of  db  with  probabllit 
density  function  l’y(.v).  bet  X be  a random  variable  in  units 
of  watts,  with  probability  density  function  K is).  Them  randi 
va  r 1 ab.l  er.  are  related  thruup;h  the  nonlinear  f ran:;  forma  t i mi : 


Y = 10  lop;  X 

1 u 


The  inverse  trans  format,  i on  is 

x - io  •v/|  0 


> 

It  Is  assumed  that  x - 0. 


Since  probalvi  1 I l y measure 
i‘e  1 at  i onr.h  i p i s valid: 


conserved  , 


the  fa  1 1 ow  i nr;  d i f fe 


l'Y(y)  dy  K (x)  ilx 


From  Equal  ion  1 and  the  differential  of  Equal,  ion  f , the  re: 
formula  for  l’y(y),  p;lveu  1"  ( x ) i : 

yio  v i o 

i’Y (y ) (o.,’.u).':»a)  10  k (10  ) 


This  r.ubreutine  i nip  1 orient,  s Equations  1 and  “’I  . The  oallinp; 
sequence  for  Lite  subroutine  is  I'Ok.'l'b  (l’X,N,l’Y)  wiiert'  N is 
the  number  of  sample  points,  in  X and  in  Y;  i'X  is  a (fxN) 
dimensioned  array  -out  a. I til  tip,  the  dens,  ity  whore 

l’X(l,J)  Fx(x)  and 

I’XC'.J)  X;  and 


nn 

(1) 


(.') 

rent  i 

( 0 

mil  i 


(»U 


'I-  1 'I 
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l’Y  Is  a (:\\N)  array  cent  a t n 1 up;  the  c>  .input  ed  values  of  the  ile 
function  In  ki i > where  l‘Y(l,J)  = I'yly)  and  PY(.’,«0  Y.  Kip.ure 
shows,  the  flowchart  for  this  subroutine.  A 1 i si  inp,  of  lie 
subroutine  is.  as  fo  lie  ws : 


sun  rout  i hk  rowrnncrx,  n,  py) 

DIMENSION  PX(2,  N)  , PY(;\  N) 

c 

C SUn  ROUT  I NK  WRITTEN  >0  NOV  loyf,  nY  7FSKTND 
C CONVERT:'  PROBAB  I I.  IT  Y DENSITY  FUNCTION  FROM  POWER  RATIO 
C TO  nn  HE:  MICRO  PASCAL 
C 

C N=  NUMBER  OF  SAMPLE  POINTS 
C PX(1  ,.!)=  P(X)  = DENSITY  IN  POWER 
C PX(2,J)=  X r POWER  RATIO 

C P Y ( 1 f J ) s P ( Y ) = DENSITY  IN  DP 
C PY(?,J)=  Y = 10  LOG  X = DB 
C 

As(10.0)«*(0.  1) 

DO  100  J=1 , N 

PY  ( 2 , J ) = 1 0.0*ALOC10(PX(?, J ) ) 

PY(1,J)=(0.2  I02BP50()2OC)iu) ) * ( A*  * P Y ( 2,  J ) ) * P X(  1 . J ) 

100  CONTINUE 
RETURN 
END 


ns.  i t y 
10 


D - IS 
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FIGURE  lb 


Flowchart  of  SUBROUTINE  POWODB. 
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4.5.3  Subroutine  DIST 


This  subroutine  compute:!  the  probability  distribution  funciton 
D (x)  from  the  input  density  function  !■’  (O.  It  is  assumed 
that  I1’  (x)  only  has  non-zero  values  forx-O.  It  is  further 
assumed  that  !•'  (x)  can  have  a uu’t  Impulse  function  6 Cxi  at 
x=0  of  magnitude  A.  The  formula  relating  DY  to  F is  given  by: 


Dx(x)  = A 6 (x ) + J Fx(C)dfJ 


Subroutine  DIST  implements  this  equation.  The  trapizoidal  method 
of  integration  is  used  to  approximate  the  integral,  with  an 
integration  step  size  depending  cn  the  sample  points  of 
Fx(x). 

The  calling  sequence  for  the  subroutine  is  DIST  ( F , N , DF , DF.LTA) 
whole  M is  the  number  of  sample  points j F is  a ( 2 x f f ) dimensioned 
array  containing  the  probability  density  function  F.  (x)  in 
row  F ( 1 , K ) and  the  corresponding  value  of  x in  row  1^ ( 2 , K ) ; 

DELTA  is  the  magnitude  of  the  unit  impulse  function;  and 
DF  is  a (d,N)  dimensioned  array  containing  the  computed  distribution 
function  D„(x)  in  row  DF(1,K)  and  the  corresponding  value  of  X 
in  row  DI-’(2,K)  . 

Figure  17  shows  the  flowchart  for  Subroutine  DIST.  The  variable 
S is  a running  sum,  initialized  to  DELTA,  which  represents  the 
integral  as  the  DO  loop  steps  through  the  values  of  X.  DELX 
is  a variable  which  is  equal  to  one-half  the  integration 
step  size  at  each  value  of  X.  A listing  of  the  subroutine 
is  as  follows: 


SUBROUTINE  DIST ( F , N, DF , DELTA ) 

DIMENSION  F (2, N) , DF(2, N) 

WRITTEN  2 DEC  1976  BY  ZESKIND 

CALCULATES  THE  DISTRIBUTION  FUNCTION  FROM  THE  INPUT 
DENSITY  FUNCTION.  DENSITY  FUNCTION  HAS  N SAMPLE  POINTS. 
TRAPIZOIDAL  INTEGRATION  . 
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Finn  re  17  Flowchart  of  Subroutine  Hist. 
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C DFLTAr  MAGNITUDE  OF  DFI.TA  FUNCTION  AT  ORIGIN  OF  F(X). 
C F( 1 , N)=  DENSITY 
C F(2,N)=  X r POWER  RATIO 
C DF ( 1 , N ) = DISTRIBUTION  FUNCTION 
C DF(  2 , N ) = X 
C 

DF ( 1 , 1 )=PFLTA 
DF(2, 1 ) = F( 2 , 1 ) 

Sr DELTA 
do  no  K=2,N 
DF(2, K)=F(2, K) 

DF.LX  = ( F(  2 , K ) -F  ( 2 , K-1  ))/2.0 
SrS  + DFLX  * ( F( 1 , K ) +F ( 1 , K- 1 )) 

DF ( 1 , K)=S 
no  CONTINUE 
RETURN 
END 


